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1.0  INTRODUCTION 


Currently,  a  number  of  different  codes  are  used  to  model  various  types 
of  propulsion  systems.  Multidimensional  models  of  solid  propellant  charges 
have  been  developed  by  Fisher  and  Graves  (1972),  Gough  (1983),  Meineke  and 
Heiser  (1989),  Groenenboom  and  Thomsen  (1989),  Fitt  et  al  (1989),  Gibelung  and 
McDonald  (1984)  and  Schmitt  (1984).  Codes  to  model  regenerative  liquid 
propellant  guns  have  been  developed  by  Steffens  et  al  (1987,  1989)  and  by 
Coffee  (1990).  The  recent  surge  of  interest  in  electrothermal -chemical  guns 
has  resulted  in  the  development  of  several  models,  including  those  by  Chen  et 
al  (1990) ,  Kashiwa  et  al  (1990) ,  Cook  et  al  (1989) ,  Winsor  and  Goldstein 
(1990),  Sinha  et  al  (1991)  and  Hsiao  et  al  (1991). 

All  these  codes  solve  sxibsets  of  a  generally  accepted  system  of  equations 
for  multiphase  flow.  They  differ  as  to  constitutive  assumptions  and  the  method 
of  solution.  Some  have  a  three-dimensional  capability.  Others  do  not. 

It  is  a  goal  of  the  US  Army  to  develop  a  Next-Generation  code  based  on 
a  flow  solver  which  is  applicable  to  all  the  types  of  gun  systems  of  current 
and  foreseeable  interest.  It  is  intended  that  the  code  have  a  three- 
dimensional  capability  and  be  able  to  address  problems  that  either  strain 
current  computational  resources  or  exceed  them.  Anticipating  that  such  a  code 
would  benefit  from  the  emerging  new  architectures  based  on  massively  parallel 
systems  of  processors,  it  is  essential  that  the  flow  solver  be  as  simple  as 
possible,  thereby  facilitating  reprogramming  to  take  advantage  of  these  new 
architectures.  Simplicity  is  also  desired  so  that  the  code  can  be  readily 
adapted  by  specialists  in  each  of  the  various  gun  system  phenomenologies ,  since 
it  is  they  who  will  propose  the  appropriate  constitutive  laws. 

Ue  have  previously  suggested  an  approach  to  the  development  of  such  a 
code  (Gough,  1991).  Ue  suggested  that  an  algorithmic  distinction  be  made 
between  the  continuous  phases  and  the  large  discrete  phases.  In  the  context 
of  a  solid  propellant  gun  the  continuous  phases  would  be  understood  to  be  the 
products  of  combustion  while  the  discrete  phases  would  be  defined  by  the 
propellant  grains.  In  our  formulation  we  did  include  the  possibility  that  the 
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continuous  phase  might  include  small  particles  or  droplets  on  the  assumption 
that  they  were  in  mechanical  equilibrium  with  the  gases.  Ve  suggested  that  a 
Continuum  Flow  Solver  (CFS)  be  developed  with  the  following  properties.  First 
and  foremost,  the  solver  had  to  be  simple  to  code  and  modify  in  order  to 
promote  portability  from  user  to  user  and  from  computer  to  computer.  Second, 
the  solver  should  ideally  be  explicit  rather  than  implicit,  since  explicit 
solvers  tend  to  be  much  simpler  to  use  and  modify,  and  interior  ballistic 
problems  rarely  require  the  sort  of  mesh  resolution  defined  by  wall  boundary 
layer  problems  which  absolutely  demand  an  implicit  treatment.  Also,  since  most 
interior  ballistic  simulations  require  wave  tracking,  integration  using  Courant 
numbers  larger  than  unity  is  not  desirable.  Third,  the  method  had  to  be 
robust.  A  shock  capturing  capability  was  desirable  in  order  to  be  able  to 
analyze  certain  classes  of  propelling  charges.  Much  more  important,  however, 
was  the  ability  to  remain  stable  in  the  presence  of  strong  porosity  gradients 
which  would  always  occur  near  the  boundaries  of  solid  propellant  charge 
increments.  Finally,  the  development  of  the  flow  solver  would  have  to 
anticipate  the  geometrical  complexities  associated  not  only  with  chamber  and 
projectile  shape  but  also  the  more  formidable  characteristics  of  the  increment 
containers  typical  of  artillery  ammunition. 

We  did  not  suggest  that  the  continuimi  solver  be  required  to  integrate  the 
motion  of  the  large  particles.  For  that  purpose  we  suggested  the  development 
of  a  second  module,  referred  to  as  a  Large  Particle  Integrator  (LPI),  to  be 
appropriately  coupled  to  the  continuum  integrator.  For  charge  designs  which 
were  liquid  based,  such  as  Regenerative  or  Bulk-Loaded  Liquid  Guns,  or  for 
certain  Electrothermal -Chemical  Guns,  the  Large  Particle  Integrator  would  not 
necessarily  be  required.  The  presence  of  droplets  could  be  represented  within 
the  continuous  phase  unless  a  separated  flow  analysis  were  justified  by  the 
availability  of  appropriate  constitutive  data. 

In  previous  studies  (Gough,  1991  and  1992),  we  have  identified  the 
principle  of  Flux-Corrected  Transport  (Boris  and  Book,  1976;  Boris  et  al,  1993) 
as  a  suitable  basis  for  the  continuum  flow  solver.  The  method  is  explicit  and 
simple  to  use.  It  has  been  shown  to  adapt  easily  to  massively  parallel  systems 


2 


(Oran  et  al,  1990).  It  has  been  shown  to  be  appropriate  to  ETC  problems  (Hsiao 
et  al,  1991)  and  preliminary  suitability  to  SP  problems  has  been  shown  by  Goug)i 
(1992). 

In  the  present  work  we  describe  two  important  steps  in-  the  development 
of  a  Next'Generation  Code,  referred  to  as  the  NGEN  Code.  We  discuss  the 
treatment  of  non-uniform  geometry  and  present  an  algorithm  for  the  Large 
Particle  Integrator  (LPI).  Non-uniform  geometry  is  treated  by  means  of  a 
structured  rectangular  mesh  in  which  only  those  cells  intercepted  by  an 
external  boundary  element  require  special  treatment.  This  approach  has  been 
developed  with  a  view  to  the  requirement  that  the  code  will  subsequently  be 
applied  to  the  simulation  of  artillery  ammunition  for  which  the  presence  of 
packaging  materials  defines  a  number  of  internal  boundary  conditions.  The 
method  described  here  is  expected  to  adapt  readily  to  that  requirement  and  to 
minimize  computational  mesh  distortion  associated  with  motion  of  the  boundaries 
as  would  occur  for  example  if  a  boundary  conformal  mesh  were  used. 

The  LPI  algorithm  is  essentially  Lagranglan,  the  motion  of  an  aggregate 
of  particles  being  tracked  explicitly.  However,  the  method  of  coupling  to  the 
continuum  flow  involves  an  attribution  of  properties,  such  as  porosity,  mass 
generation  and  so  forth,  whose  spatial  variation  can  be  explicitly  controlled 
near  the  boundaries  of  each  charge  increment,  thereby  defining  sufficiently 
smooth  distributions  to  maintain  stability  of  the  continuum  solver. 

In  Chapter  2.0  we  discuss  the  governing  equations  presently  encoded.  The 
method  of  solution  is  presented  in  Chapter  3.0.  In  Chapter  4.0  we  illustrate 
the  code  capabilities  by  reference  to  two  data  bases,  one  for  a  155mm  howitzer 
and  one  for  a  120mm  tank  gun.  The  discussion  of  Chapter  3.0  is  supported  by 
Appendix  A  which  presents  a  characteristic  analysis  of  the  equations  of  two- 
dimensional  inviscid  flow.  Appendix  B  presents  a  listing  of  the  current  NGEN 
input  files . 
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2.0  THE  GOVERNING  EQUATIONS 

Ue  first  provide  a  general  statement  of  the  balance  equations  for  a 
multiphase  mixture  In  Section  2.1.  In  Section  2.2  we  present  the  reduced  set 
idilch  Is  encoded  Into  the  present  version  of  the  NGEN  Code.  In  Section  2.3  we 
discuss  the  current  set  of  constitutive  laws.  In  Section  2.4  we  note  the 
boundary  conditions  and  the  equation  of  motion  of  the  projectile. 

2.1  General  Statement  of  the  Balance  Equations 

The  various  systems  of  equations  which  have  been  considered  in  current 
models  are  all  subsets  of  the  macroscopic  balance  equations  for  a  multi- 
component  mixture  as  discussed  elsewhere  (Gough,  1992). 

The  mixture  may  be  viewed  as  consisting  of  a  multicomponent  fluid, 
referred  to  as  the  continuous  phase,  and  an  aggregate  of  droplets  or  solid 
particles,  referred  to  as  the  discrete  phase.  The  continuous  phase  is 
understood  to  be  a  multi- component  mixture  of  gases  and  droplets  or  particles 
which  are  small  enough  to  ensure  local  mechanical  equilibrium.  The  gases  are 
always  in  local  thermal  equilibrium  while  the  droplets  and  particles  are  not 
necessarily  so.  The  continuous  phase  is  characterized  by  single  local  values 
of  density,  p\  velocity  vector,  u;  temperature,  T;  pressure,  p;  shear  stress 
tensor,  r;  and  internal  energy,  e.  It  Is  assumed  to  comprise  N^  species  each 

characterized  by  local  values  of  mass  fraction  ,  1-  1 . ,Nj.  Moreover, 

the  velocity  u  is  understood  to  be  the  barycentric  or  mass  wel^ted  average  of 
the  velocities  of  each  of  the  components  (Williams,  1965).  Each  component  is 
characterized  by  a  diffusion  velocity  v^  relative  to  the  barycentric  value  u. 

The  term  discrete  phase  is  understood  to  refer  to  an  aggregate  of 
particles  or  droplets.  If  a  solid  propellant  charge  is  being  modeled,  each 
type  of  propellant  will  constitute  a  component  of  the  discrete  phase.  Other 
components  may  be  present  if  the  decomposition  of  the  propellant  or  an  igniter 
element  yields  intermediate  combustion  products  in  particulate  or  droplet  form. 
Still  other  components  may  need  to  be  considered  if  rupture  of  a  container  is 
to  be  modeled  or  if  wear-reducing  additives  like  talc  are  present  and  their 
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dispersal  pattern  Is  to  be  calculated.  In  the  case  of  the  RLPG  the  aggregate 
may  consist  of  a  spray  created  by  the  breakup  of  the  Injected  Jet  of  liquid 
propellant.  Similarly,  In  the  ETC  the  aggregate  may  consist  of  droplets 
created  by  the  Helmholtz  instability  on  the  boundary  of  the  Taylor  cavity. 
Because  of  the  generality  Implicit  In  the  representation  of  the  continuous 
phase,  the  foregoing  components  of  the  discrete  phase,  except  for  the  solid 
propellant  grains,  may  be  optionally  Included  In  the  continuous  phase  and 
modeled  according  to  a  homogeneous  mixture  representation.  Alternatively,  they 
may  be  modeled  Independently  when  the  situation  so  warrants. 


Ve  assume  that  the  discrete  phase  consists  of  a  total  of  components. 
Each  component  is  characterized  by  density,  stress  tensor,  a ;  velocity 

vector,  Uj^;  temperature,  number  density,  i:\j^ ;  and  morphological  data  to 
characterize  the  volxame,  ,  and  surface  area,  of  each  particle  or  droplet. 
The  temperature  T^^  may  be  either  a  surface  or  a  bulk  temperature  depending  on 
the  nature  of  the  model  assumptions  for  the  problem  In  question. 


In  general  It  Is  necessary  to  consider  a  macroscopic  formulation  of  the 
balance  equations  due  to  the  presence  of  the  discrete  phases.  The  macroscopic 
formulation  is  such  that  it  reduces  to  standard  continuum  equations  in  the 
context  of  single-phase  flow.  Given  a  microflow  property  we  use  <V^>, 
and  <^^>8  to  respectively  denote  a  bulk  average,  a  mass-weighted  (Favre) 
average  and  an  interphase  surface  average. 


We  use  a  to  denote  the  porosity,  or  the  fraction  of  a  unit  macroscopic 
volume  occupied  by  the  continuous  phase.  Similarly,  we  use  to  denote  the 
volume  fraction  of  the  i-th  discrete  phase.  Evidently, 


Q 


Hd 


1  -  E 


2.1.1 


Ue  may  now  state  the  balance  equations  for  the  continuous  phase  in  the 
following  forms  which  differ  from  those  presented  earlier  (Gough,  1992)  only 
in  respect  to  the  neglect  of  certain  correlation  terms  which  result  from  the 
formal  macroscopic  averaging  process. 
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Ve  have  the  balance  of  mass 


Hd 


—  la<p>]  *  V»[a/)<u>^]  -  ^ 

Bt  i.i 


2.1.2 


where  m^^  is  the  rate  of  decomposition  per  unit  volume  of  the  i-th  discrete 
phase . 


Each  of  the  j  contponents  of  the  continuous  phase  satisfies  a  balance  of 
mass  equation  in  the  following  form, 

.±[a<p><yj>^]  +  V*[a<p><Yj>^{<u>p  ♦<Vj>^)] 


-  Q<«j>  ♦  ^  md^<Yd^j>, 
i-1 


Here  <Wj>  is  the  average  rate  of  production  per  unit  volume  of  species  j  by 
chemical  reactions  and  is  the  average  mass  fraction  of  species  j 

produced  by  the  decomposition  of  the  i-th  discrete  phase  including  the  effect 
of  the  surface  reaction.  The  macroscopic  balance  of  momentum  for  the 
continuous  phase  takes  the  form 

[o<p><u>p]  +  V*  [o<p><u>^<u>p] 

5t 

Rd  Rd 

-  -  aV<p>  +  oV‘<r>  ♦  md^<Udj>  -  nd^Sd^fdj  2.1.4 

1-1 

Here  f^^  represents  the  interphase  drag  due  to  the  i-th  discrete  phase. 
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The  energy  equation  for  the  continuous  phase  takes  the  form 
^  |q<P>  j^<e>p  +  1<U>*  j  j.  V*  |q  <,><u>.  |<e>,  . 


V*  Q<u>p  •  <a>  -  V*  Q<q>  -  <p>  ndiSdi<qdi>. 

Bt  1.1 


••d 


2.1.5 


Hd  "d 

E  ndiSdifdi'<«di>p  *  E 


1«1 


<ed^>  ♦  52^.  ♦ 


<«di>^ 


<Pdi> 


We  note  the  phase  Interaction  terms  on  the  right  hand  side  of  Equation  2.1.5. 
The  first  of  these  introduces  the  heat  transfer  due  to  conduction  and  radiation 
per  unit  surface  area  The  next  term  represents  the  work  done  by  the 

interphase  drag.  The  third  term  represents  the  heat  added  due  to  decomposition 
of  the  discrete  phases. 

The  i-th  discrete  phase  is  found  to  be  governed  by  a  macroscopic  mass 
balance  analogous  to  that  for  the  continuous  phase,  namely 

•^[odi<Pdi>]  +  V*[“di<Pdi><'idi>p]  =  “  “di  •  2.1.6 


The  i-th  discrete  phase  is  governed  by  a  macroscopic  balance  equation  in 
the  following  form 


-^[“di<Pdi><'idi>pl  ^•[“di<Pdi><«di><«di>] 

“  -adi^<P>  - 


2.1.7 


We  note  on  the  right  hand  side  of  Equation  2.1.7  the  formal  presence  of  a 
stress  term  which  reflects  the  difference  between  the  average  stress  in 

the  i-th  discrete  phase  and  the  average  ambient  pressure  in  the  continuous 
phase.  This  is  interpreted  as  reflecting  interactions  between  droplets  or 
particles. 
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2.2  Balance  Equations  for  Current  NGEN  Code 


In  the  current  version  of  the  code  diffusive  processes  are  not  consider¬ 
ed.  We  consider  the  mixture  to  consist  of  a  continuous  phase  -  the  combustion 
products  -  with  constant  thermodynamic  properties,  and  a  single  discrete 
aggregate  -  the  propelling  charge.  In  order  to  provide  a  vehicle  for  the 
representation  of  an  igniter  we  Include  the  Influence  of  a  local  source  term. 
Ve  consider  the  continuous  phase  first.  Since  there  is  no  composition 
dependence,  Equation  2.1.3  is  not  required.  Moreover,  we  drop  the  explicit 
representation  of  the  averages.  The  resulting  balance  equations  for  the 
continuous  phase  are  therefore 


gap 

at 


dz 


apu  -  — 


r  ar 


2.2.1 


aapu  a  la  ao*  ^ 

- —  “  -  — _  apuu  -  Jz  JL  ropuv  -  ag_  ^  ♦mu  -  f 

at  az  r  ar  az  * 


2.2.2 


apvu  -  1.  rapvv  -  ag  ^2.  ♦  mv  -  f 

at  az  r  ar  ar 


2.2.3 
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ao 
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U  “  +  v; 


2gc 


+  ®i«e 


i8®ig 


.  2.2.4 


where  E  ■  p  I  e  ♦  JLLH  I  and  g^  is  used  to  reconcile  units  and  we  use  a  subscript 

I  2g„  J 

p  to  denote  a  property  of  the  solid  propellant  which  constitutes  the  discrete 
phase.  Also  m  and  qp  are  respectively  the  rates  of  mass  and  heat  transfer  per 
unit  volume  while  f,  and  f,  are  the  components  of  the  interphase  drag.  Ve  use 
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to  denote  the  rate  of  addition  of  a  source  term  per  unit  volume  and  e^g 
represents  the  chemical  energy  of  the  source. 

-As  we  will  discuss  further,  in  Chapter  3.0,  these  Equations  are 
Integrated  using  the  one -dimensional  solver  LCPFCT  which  assumes  that  the 
system  is  timesplit  Into  an  axial  set 

l££  -  -  -Lapu  ♦  m  ♦  m^  ,  2.2.5 

at  dz 

—  opuu  -  ag^  .££  ♦  mu  -  f,  ,  2.2.6 

at  dz  dz 


and  a  radial  set 


8ap  . 

-  1  J.  rav  , 

2.2.9 
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r  dr 

dapu 

1  fl 

*  -  _  rapuv  , 

2.2.10 

at 

r  ar 

dapv 

»  -  _  rapvv  -  ag 

^  +  iv  -  f^ 

1 

2.2.11 

at 

r  ar 

ar 
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-  ■  - - ravE  -  —  —  rapv  .  <>  o  no 

at  r  ar  r  ar 

The  balance  equations  for  the  discrete  phase  •  the  aggregate  of  granular 
propellant  •  may  be  restated  in  the  present  context  as  follows.  Equation  2.1.6 
becomes 


-1(1  -o)  ♦V'fd  -a)Up]  •  -  JL  .  2.2.13 

3b  Pp 

Equation  2.1.7  may  be  simplified  with  the  help  of  2.1.6  and  the  assumption  that 
is  Isotropic  to  yield 

Du 

p  (l-o) - P  -  -g„(l  -  o)Vp  -  g„7(7  +  f  ,  2.2.14 

where  D/Dtp  is  the  convective  derivative  along  the  propellant  streamline.  We 
do  not  express  2.2.13  and  2.2.14  in  cylindrical  coordinates  because,  as  we 
discuss  in  Chapter  3.0,  they  are  solved  in  a  Lagranglan  formulation  for  which 
2.2.14  is  the  natural  expression.  We  note  that  o  in  Equation  2.2.14 
corresponds  to  in  Equation  2.1.7. 

2.3  Constitutive  Laws 

Closure  of  the  balance  equations  requires  equations  of  state  for  the 
continuous  phase  and  the  discrete  phase  as  well  as  correlations  to  describe  the 
interphase  transfer  processes  -  heat  transfer,  drag  and  combustion  -  together 
with  laws  to  determine  the  morphology  of  the  discrete  phase  -  the  surface  area 
and  volume  of  the  particles.  In  general,  the  laws  chosen  here  are  a  subset  of 
those  used  in  current  versions  of  the  TDNOVA  Code  (Gough,  1983). 
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2.3.1  Equation  of  State  of  Continuous  Phase 


It  is  asstoned  that  the  continuous  phase  obeys  the  covolune  equation  of 

state 

e  -  c„T  -  ,  2. 3. 1.1 

(7-l)p 

where  b  is  the  covolume,  7  is  the  ratio  of  specific  heats  and  Cy  is  the 
specific  heat  at  constant  volxune. 


The  molecular  weight  and  the  ratio  of  specific  heats  are  assumed  to  be 
constant  and  are  given  values  appropriate  to  the  fully  reacted  propellant. 

2.3.2  Equation  of  State  for  Discrete  Phase  (Intergranular  Stress  Law) 


The  intergranular  stress  is  taken  to  depend  on  porosity  and  also  on  the 
direction  of  loading.  Ue  embed  the  constitutive  law  into  the  formula  for  the 
rate  of  propagation  of  Intergranular  disturbances 


a(c) 


go  do 
Pp  <kt 


2. 3. 2.1 


Ue  may  recast  2. 3. 2.1  into  a  form  more  suitable  for  numerical  inte¬ 
gration,  namely 


--Pp±!-^ 

go  Dtp 


2. 3. 2. 2 


In  order  to  formulate  the  functional  behavior  of  a(a)  we  introduce  ,  the 
settling  porosity  of  the  bed,  and  values  of  a(a)  equal  to  a^  and  a2  which 
respectively  correspond  to  loading  at  and  to  unloading/reloading.  The 
nominal  loading  curve,  corresponding  to  monotonlc  compaction  of  the  bed  from 
Qg  to  a  smaller  value  of  the  porosity  a  is  given  by 


O  •  On<»(o) 


2. 3. 2. 3 
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The  functional  dependence  of  a(a)  may  now  be  stated  as: 


a(o) 


aiQ^/a  if  o  ^  0.  o  -  o  ^ 

a^  if  0  ^  <7<  OnoB.  a  ^  Oo 

or  if  d  >  0,  ff  »  OnoBo  “  ^  ®o 


2. 3. 2. 4 


0  if  a  “  0  and  d  >  0  or  if  o  >  Oq 


where  we  understand  d  to  mean  Da/Dtp. 

^ . 3.. ? _ Discrete  Phase  Morphology  (Fomi  Functions) 

It  is  assumed,  in  the  present  study,  that  the  propellant  grains  are 
multi-perforated  cylinders  having  initial  length  ,  external  diameter  and 
perforation  diameter  d^ .  Until  such  time  as  slivering  occurs,  that  is  to  say 
the  time  at  which  the  regressing  perforation  surfaces  intersect,  the  surface 
area  and  volume  are  given  by 

Sp  .  k(L„  -  2d)l(Do  -  2d)  +  N(d„  +  2d)] 

2. 3. 3.1 

+  ir/2[(D„  -  2d)2  -  N(d„  *  2d)2] 

Vp  -  ir(L^  -  2d)((D«  -  2d)2  -  N(d„  ♦  2d)2]/4  ,  2. 3. 3. 2 

where  N  is  the  number  of  perforations  and  d  is  the  total  linear  surface 
regression,  assumed  uniform  over  all  the  surfaces  of  a  grain. 
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Once  slivering  occurs,  the  form  functions  become  rather  complicated  for 
N  >  1.  Formulas  for  the  form  functions  following  the  slivering  of  seven- 
perforation  grains  may  be  found  in  Krier  et  al  (1973).  The  present  version  of 
the  code  supports  single-,  seven-  and  nineteen-  perforations  grains  in  the 
propelling  charge. 

2 .3. A  Interohase  Drag 

We  express  the  interphase  drag  in  a  granular  bed  in  the  following  form 
(Gough,  1983) 

f.  -  W"  pju -Up|(u  -Up)f,  2. 3. 4.1 

l“eJ 


Equation  2. 3. 4.1  refers  to  the  exterior  voidage,  a,,  and  the  effective 
diameter,  ,  based  on  the  exterior  volume  and  surface  area  computed  as  thoug)i 
the  perforations  were  not  present. 


f. 


ci 


»,R6 


max 


Ci, 


a.RG 


1  -  C.  ".o 


1  -  c,  c, 

®o  • 


0.45 


min 


if  «,  ^  c,^ 

2. 3. 4. 2 

if  <  c  ^  1  , 


C 


r 

0.85  if  the  grain  is  perforated  and  unignited 
1.0  otherwise  , 


2. 3. 4. 3 


f 


min 


{0.3  for  spheres 
0.75  for  cylinders 


2. 3. 4. 4 


t 


•  ,RG 


2.50  Re;'’  ®®^A2  i7 


2. 3. 4. 5 


2. 3.4.6 


13 


1 


for  spheres 


2. 3. 4. 7 


A 


[(1/2  ♦  L/D>/(3L/2D>2/3 


for  cylinders 


We  have  used  as  the  value  of  o,  in  the  settled  condition  of  the  bed, 

L  as  the  length  of  a  cylindrical  grain  and  D  as  its  diameter,  and  and  as 
the  density  and  viscosity  of  the  gas  at  the  film  temperature. 

2.3.5  Intemhase  Heat  Transfer 

The  interphase  heat  transfer,  in  both  the  propelling  charge  and  the 
centercore  ignition  charge  is  assumed  to  be  governed  by  the  empirical 
correlation  of  Gelperin  and  Einstein  (1971) .  We  express  the  heat  transfer  in 
the  fora 

NUp  *0.4  Pr^^®Rep*^®  ,  2. 3. 5.1 

where 


NUp  -  hDp/kf 

Rep  -  p,|u  -Up|Dp//i£ 
h  =  q/(T-Tp) 

where  h  is  the  film  coefficient  and  q  is  the  heat  transfer  per  unit  surface 
area.  The  subscript  f  denotes  an  evaluation  of  properties  at  the  film 
temperature  (T  +  Tp)/2  where  T  and  Tp  are  respectively  the  continuous  phase  bulk 
average  temperature  and  the  particle  surface  average  temperature.  The 
viscosity  is  taken  to  have  a  Sutherland- type  dependence  on  temperature, 

M  -  0.134064  2. 3. 5. 2 

T  4  110 
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The  thenaal  conductivity  follows  from  the  Prandtl  number  which  is  assumed  to 
satisfy 


Pr  -  1st  •  — fLZ_  2. 3. 5. 3 

k  97-5 

The  heat  transfer  per  unit  volume  qp  is  related  to  q  according  to 

qp  ■  (1  -  o)  is  q  2. 3. 5. 4 


2.3.6  Solid  Phase  Surface  Temperature 


Assuming  that  ignition  is  an  essentially  uniform  event  with  respect  to 
the  surface  of  each  grain  of  either  the  propelling  charge  or  the  centercore 
ignition  charge  and  supposing  that  the  temperature  distribution  within  the 
solid  phase  can  be  captured  by  a  cubic  profile,  leads  to  the  following 
expression  for  the  surface  temperature 


*Po 


2  hH 


2. 3. 6.1 


where  Tp^  is  the  initial  surface  temperature  and  H  satisfies 

*  o  q  2  3  6  2 

Dtp  "  ■ 

2.3.7  Ignition  and  Combustion 

Ignition  is  assumed  to  occur  when  the  surface  temperature  exceeds  a  pre 
determined  value.  The  rate  of  surface  regression  is  given  by 
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Bj  ♦  BjP® 


2. 3. 7.1 


Dd 


It  should  be  noted  that  only  one  of  2. 3. 6. 2  and  2. 3. 7.1  has  to  be  solved  at 
each  point  according  as  the  temperature  is  less  than  or  equal  to  the  ignition 
temperature . 


The  mass  transfer  per  unit  volume,  m,  is  related  to  the  regression  rate 
d  according  to 


m  -  (1  -  a)^d 


2. 3. 7. 2 


Ui _ B_ouiidarv_Condltions  and  Projectile  Motion 

For  the  continuous  phase  we  impose  slip  boundary  conditions  at  all 
bounding  surfaces.  For  the  discrete  phase  we  admit  the  possibility  of 
separation  from  an  external  boundary  so  that  the  appropriate  condition  is  one 
of  non-penetration. 

The  projectile  is  assumed  to  move  as  a  rigid  body  subject  to  the  total 
force  on  its  base  and  the  surface  of  an  afterbody,  if  present,  and  the 
resistance  due  to  interactions  with  the  gun  tube.  The  latter  may  be  expressed 
in  a  variety  of  forms.  Appendix  B  may  be  consulted. 


16 


3.0  METHOD  OF  SOLUTION 


As  we  have  already  noted,  we  have  selected  the  Flux- Corrected  Transport 
Algorithm  LCPFCT  (Boris  et  al,  1993)  as  the  basis  for  the  continuum  flow 
solver.  In  Section  3.1  we  briefly  siammarlze  the  algorithm  and  note  Its 
application  at  boundary  points  and  its  treatment  of  multidimensional  flows. 
In  Section  3.2  we  turn  to  the  problem  of  non-uniform  boundary  geometry.  We 
discuss  the  use  of  the  Virtual  Cell  Principle  (Landsberg  et  al,  1993)  as  a 
means  of  accommodating  arbitrary  boundary  geometries  within  a  uniform, 
structured,  rectangular  mesh.  We  compare  the  Virtual  Cell  approach  with  the 
method  ultimately  adopted  in  the  NGEN  Code,  namely  one  in  which  the  actual 
boundary  cell  geometry  is  used  and  the  solution  stabilized  by  a  characteristic 
analysis  of  the  boundary  values.  The  comparison  is  made  in  the  context  of  one- 
dimensional  flow.  In  Section  3.3  we  discuss  the  implementation  of  the 
characteristic  based  boundary  analysis  for  two-dimensional  flow.  In  Section 
3.4  we  discuss  the  Large  Particle  Integrator. 


3.1  Suinmarv  of  LCPFCT  Algorithm 


LCPFCT  is  a  one -dimensional  solver  for  a  canonical  balance  equation  in 
the  form 


^  -  -L-  ±(rk-ipv)  -  -L-  J_(rk-iDi) 


3.1.1 


where  k  -  1,  2  or  3  for  planar,  cylindrical  and  spherical  flows  respectively. 
Here  p  is  the  transported  variable  and  ,  Dj  and  Dj  are  referred  to  as  source 
terms . 


The  computational  domain  is  decomposed  into  N  cells .  The  state  variable 
Pi  is  considered  to  apply  to  the  center  of  the  i-th  cell.  The  cell  has  volume 
Ai  and  is  bounded  by  surfaces  whose  areas  are  A^.^/a  which  are 
presumed  to  be  orthogonal  to  the  fluid  streamline.  Fluid  properties  on  the 
cell  boundaries  are  determined  by  averaging  with  the  values  in  adjacent  cells. 
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The  fluid  velocity  at  the  right  hand  boundary  of  cell  i  is  while  the 
velocity  of  the  boundary  is  .  We  define  Av^^^^2  ”''^1*1/2  "  ''^i»i/2* 

Equation  3.1.1  is  integrated  according  to  a  finite  volvime  formulation  via 
several  steps  in  which  strong  diffusion  is  first  introduced  and  then 
subsequently  canceled  by  antidiffusion  to  the  maximum  extent  consistent  with 
the  minimal  requirements  of  computational  stability  and  the  condition  that  the 
antidiffusion  not  introduce  new  maxima  or  minima  in  the  updated  distribution 
of  p.  Properties  at  the  beginning  of  the  timestep  are  denoted  by  superscript 
o  while  updated  properties  are  denoted  by  superscript  n. 


We  first  discuss  the  scheme  for  non  boundary  points .  The  treatment  of 
boundary  points  is  taken  up  subsequently.  First  p  is  transported  over  timestep 
At  to  define 


AiPi  =  AiP^  -  AtPi^j/2  Ai^j/2  AVi^j/2  +  AtPi_j/2  A^_j/2  AVi_j^/2  .  3.1.2 

Next,  the  effect  of  the  source  terms  is  added  to  define  p^  according  to 


Aj  p  =  A^p^  +  _  A^^j^2  (  *  ^i,i  )  “  —  (  ^i,i  ^1,1-1 ) 

^i-1/2  )  {  ®2,i*l  "  ®2,l-l) 

4 


3.1.3 


Except  for  the  influence  of  the  change  in  cell  volume  from  A°  to  a”  ,  Equations 
3.1.2  and  3.1.3  represent  the  discretization  of  Equation  3.1.1  according  to  a 
finite  volume  formulation.  The  effect  of  the  change  in  cell  volume  is  combined 
with  a  strongly  diffusive  process  according  to 


A"?, 


.0  T 
A.p.  *  V 


i-l/2 


1/2  ^i-1/2 


(  '  P\-x  ) 


3.1.4 


where  A^^j^2 


and  1/  is  a  coefficient  of  numerical  diffusion. 
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Finally,  as  discussed  by  Boris  et  al  (1993),  flux*corrected  antldlffuslve 
fluxes  are  computed  and  the  final  updated  value  of  p  Is  given  by 


Boundary  conditions  are  supplied  at  ^v*uz  ’  formulated 

In  conjunction  with  guard  cells  located  outside  the  boundary  so  that  the 
Integration  of  boundary  cells  can  be  performed  using  essentially  the  same  code 
as  that  which  Is  used  for  the  Interior  cells.  In  general  LCPFCT  expects  the 
user  to  specify  the  boundary  motion  and  sets  the  left  hand  boundary  velocity 
difference 

^^1/2  “^f/2  -  -  rr/2)/^t  .  3.1.6 

and  AVj,^j^2  Is  defined  analogously.  For  nonperiodic  boundary  conditions  the 
guard  cell  value  pg  Is  permitted  to  take  the  general  form 

Pq  *  ^  ®2  •  3.1.7 

where  B^,  B2  are  set  by  the  user.  The  right  hand  guard  cell  value  p^  is 

defined  analogously.  The  report  by  Boris  et  al  (1993)  may  be  consulted  for  a 
general  discussion  of  the  choices  of  B^  and  Bj  for  a  variety  of  boundary  types. 

As  discussed  by  Boris  et  al  (1993)  LCPFCT  may  be  made  second  order 
accurate  in  time  by  using  a  two-step  scheme,  first  integrating  from  t  to  t  + 
At/2  and  then  from  t  to  t  -t-  At  using  the  intermediate  values  to  define  the 
geometric  terms  and  the  source  terms.  Further,  LCPFCT  Is  applied  to 
multidimensional  situations  by  timesplitting,  integrating  first  in  one 
direction  and  then  another  with  the  appropriate  subsets  of  the  multi¬ 
dimensional  equations. 

Then  referring  to  the  tlmespllt  system  of  equations  for  the  continuous 
phase.  Equations  2.2.5  •  2.2.12,  It  Is  easy  to  see  that  each  of  these  Is  of  the 
form  of  Equation  3.1.1  with  the  proper  choices  of  p,  ,  D2 ,  D3  and  k. 

However,  while  a  natural  set  of  choices  for  the  fluxes  would  be  ap,  apu,  apv 
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and  qE,  we  have  elected  to  use  p,  pu,  pv  and  E,  folding  the  influence  of  the 
porosity  a  into  the  metric  quantities  This  selection  yields 
generally  smooth  thermodynamic  state  variables  and  minimizes  the  influence  of 
the  diffusive  step  defined  by  Equation  3.1,4. 

As  for  the  selection  of  guard  cell  values,  Equation  3.1.7,  we  are 
concerned  here  only  with  impermeable  boundaries  with  slip.  We  note  that  when 
the  boxindary  is  impermeable  the  guard  cell  values  only  affect  the  diffusion 
terms.  We  therefore  use  the  simplest  choices  of  and  B2 ,  namely  those  which 
provide  even  reflection  of  the  thermodynamic  properties  and  the  tangential 
velocity  and  odd  reflection  of  the  normal  velocity  relative  to  the  boundary. 

3.2  Treatment  of  Non>Uniforni  Geometirv 

It  is  immediately  obvious  how  to  use  LCPFCT  to  integrate  multidimensional 
flows  on  domains  whose  boundaries  are  rectangular  and  conformal  with  the 
coordinate  surfaces  of  a  Cartesian  grid.  Here  we  consider  the  application  of 
the  method  to  non-uniform  geometry. 

Much  progress  has  been  made  in  recent  years  with  regard  to  the  analysis 
of  flows  with  arbitrary  complex  three-dimensional  boundary  conditions.  Methods 
based  on  unstructured  meshes  have  been  shown  to  be  especially  powerful  for 
addressing  complex  geometries,  see  for  example  Lohner  (1989),  but  do  not  adapt 
well  to  massively  parallel  computing  systems.  The  popular  mapping  algorithms 
of  the  type  discussed  by  Thompson  et  al  (1985)  achieve  block  structured  arrays, 
but  are  not  as  flexible  as  the  unstructured  mesh  algorithms.  Moreover,  these 
methods  result  in  non-uniform  mesh  lines  which  introduce  computational  errors 
through  the  need  to  resolve  accurately  the  metric  derivatives  associated  with 
the  transfoxnnation.  The  non-uniformity  can  also  result  in  regions  where  the 
mesh  is  unnecessarily  crowded  with  the  result  that  the  timestep  is  reduced  in 
order  to  satisfy  the  stability  condition  for  an  explicit  integration  scheme. 

Recently,  Landsberg  et  al  (1993)  have  suggested  an  alternative  to  these 
boundary  fitted  mesh  algorithms.  They  suggest  the  use  of  a  simple  rectangular 
mesh  which  covers  a  superset  of  the  computational  domain  and  extends  past  the 
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physical  boundaries  In  various  locations.  Cells  are  either  open,  closed  or 
partially  obstructed  according  as  they  lie  Inside  the  boundaries,  outside  the 
boundaries  or  over  the  boundaries.  Only  the  partially  obstructed  cells  require 
special  treatment.  According  to  Landsberg  et  al,  accurate  solutions  can  be 
obtained  provided  that  the  partially  obstructed  cells  are  characterized  by 
correct  values  of  volume,  face  area  and  the  normal  vector  to  the  physical 
boundary.  The  normal  vector  Is  used  to  determine  the  flux-coupling  which 
results  when  the  balance  equations  are  tlmespllt  as  discussed  previously. 
Since  the  cell  faces  are  formally  located  the  same  distance  apart  as  those 
which  are  open,  the  stability  condition  does  not  demand  a  reduced  tlmestep. 
Some  distortion  of  the  flow  Is  expected  near  the  boundary  since  the 
Interpolation  of  flow  properties  Is  not  consistent  with  the  fraction  of  the 
cell  which  Is  actually  Inside  the  computational  domain.  As  discussed  by 
Landsberg  et  al,  sufficiently  accurate  characterizations  of  cell  volume, 
surface  area  and  average  normal  vector  can  be  determined  by  subdividing  the 
boundary  cell  Into  "virtual  cells"  whose  purpose  Is  simply  to  establish  these 
metric  data  and  for  which  state  variables  are  never  determined. 

In  the  present  application,  the  metric  data  can  be  determined  by  simple 
analytical  formulas,  so  that  cell  subdivision  Is  not  required.  Nevertheless, 
we  will  refer  to  the  partially  obstructed  cells  whose  boundaries  lie  outside 
the  computational  domain  as  virtual  cells. 

The  Virtual  Cell  Principle  has  the  advantages  that  It  allows  easy 
parallelization  of  the  algorithm,  complete  flexibility  with  regard  to  boundary 
geometry,  minimal  additional  computation  at  the  boundaries,  and  freedom  from 
the  expense  and  possible  inaccuracies  associated  with  the  computation  of  non- 
uniform  mesh  transformation  derivatives. 

Accordingly,  we  examine  the  applicability  of  the  Virtual  Cell  Principle 
to  the  NGEN  Code,  using  as  a  baseline  the  analytical  solution  of  Love  and 
Pidduck  (1921)  of  the  well-known  Lagrange  problem. 

At  the  same  time  we  evaluate  an  alternative  procedure  for  the  boundary 
cells.  Ve  again  take  the  mesh  to  be  orthogonal  and  rectangular.  Cells  are 
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again  open,  closed  or  partially  obstructed.  But  for  the  partially  obstructed 
cells  we  analyze  the  physical  element  defined  by  the  fraction  of  the  cell  which 
is  actually  inside  the  computational  domain.  Thus  whereas  the  virtual  cell 
analysis  only  corrects  Ak  and  partially  distributed  cells,  the 
physical  cell  analysis  also  corrects  location  of  the  boundary.  This 
raises  the  problem  of  violation  of  the  stability  condition  and  we  circumvent 
this  by  means  of  an  implicit  flux  between  the  boundary  cell  and  its  Internal 
nel^bor  which  forces  the  updated  states  to  satisfy  a  characteristic  condition 
of  compatibility  based  on  states  N  •  1,  N  and  the  physical  boundary  condition. 
Ve  refer  to  this  approach  as  the  characteristic  based  method.  The  two 
procedures  are  compared  in  Figure  3.1. 


VIRTUAL  CELL  PHYSICAL  CELL 


Volume:  A„  =  f 


Ah  -  f 


Area: 


Base  Pressure:  Pb  “  ( f  “  -)  ( Pn  "  -  i)  ^  -  i ) 

2  1  +  f 

Figure  3.1.  Comparison  of  Virtual  Cell  and  Physical  Cell  Analysis  of 
Boundary  Element  for  One •Dimensional  Problem. 


In  the  characteristic  based  method  the  value  of  Pb  as  computed  from  the 
interpolation  formula  based  on  the  LCPFCT  update  Is  compared  with  the  value 
based  on  a  characteristic  analysis  of  data  at  N  •  1  and  the  boundary,  where  the 
velocity  is  known.  Cells  N  -  1  and  N  then  exchange  a  quantity  of  mass  which 
brings  the  two  values  into  agreement. 
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For  this  simple  one -dimensional  flow  the  appropriate  characteristic 
condition  may  be  stated  as  (Courant  and  Hilbert,  1953} 
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which  is  appropriate  at  the  right  hand  boundary,  the  base  of  the  projectile. 
At  that  location  du/dt  Is  known  from  the  motion  of  the  projectile.  The  spatial 
derivatives  dp/dz  and  du/dz  can  be  discretized  using  the  old  data  at  cell 
N  -  1  and  the  old  boundary  values,  the  boundary  value  of  pressure  being 
determined  by  extrapolation  as  shown  In  Figure  3.1.  Then  3.2.1  allows  the 
determination  of  the  characteristic  based  value  of  pressure  at  the  new  time. 


The  benchmark  problem  used  for  the  evaluation  of  the  two  methods  entails 
motion  of  a  projectile  through  a  cylindrical  tube  due  to  the  expansion  of  a 
region  of  pressurized  gas .  An  analytical  solution  was  obtained  by  Love  and 
Pidduck  for  specific  conditions  (1921)  and  has  been  used  as  a  benchmark  in  the 
evaluation  of  other  interior  ballistic  codes  (Schmitt,  1981;  Robbins,  1983). 
The  tube  diameter  is  15  cm.  The  Initial  length  of  the  gas  column  is  169.8  cm 
and  the  projectile  travel  is  600  cm.  The  projectile  mass  is  50  kg.  The 
initial  pressure  and  temperature  of  the  gas  are  621.09  MPa  and  2666.8  K.  The 
gas  has  molecular  weight  23.8  g/gmol,  covolume  1.0  cmVg  and  ratio  of  specific 
heats  equal  to  1.22. 

To  examine  the  influence  of  the  fraction  of  a  cell  outside  the  boundary 
we  construct  an  accordion  mesh  which  expands  with  the  projectile  and  allows  a 
fraction  OUTFR  of  the  last  cell  to  be  in  front  of  the  projectile.  Here  OUTFR 
-  1  -  C  where  f  is  as  shown  in  Figure  3.1.  Ve  obtained  solutions  with  15,  30 
and  60  cells  and  various  values  of  OUTFR  ranging  from  0  (cell  completely 
inside)  to  0.99  (cell  almost  completely  outside).  Guard  cell  values.  Equation 
3.2.7,  were  chosen  to  give  symmetry  in  the  density  and  energy  and  antisymmetry 
in  the  velocity. 

It  was  found  that  both  methods  provided  results  at  muzzle  exit  which  were 
in  good  agreement  with  the  analytical  solution,  with  the  physical  element 


23 


analysis  being  somewhat  better,  especially  for  large  values  of  OUTFR,  as  might 
be  expected.  However,  early  In  the  solution,  the  virtual  cell  analysis 
produced  large  wiggles  and  errors  In  the  boundary  values.  Following  some 
review,  -It  was  found  that  by  replacing  the  velocity  V],  by  an  interpolation 
between  N  •  1  and  the  boundary,  more  stable  results  could  be  obtained.  This 
measure  could  have  been  made  contingent  on  OUTFR,  but  for  simplicity  was 
adopted  uniformly.  It  was  also  applied  to  the  characteristic  based  method 
%/here  It  had  a  slightly  beneficial  effect  In  that  It  eliminated  a  pressure 
wiggle  of  magnitude  about  IZ  near  the  forward  boundary. 

The  results  at  muzzle  exit  obtained  with  the  two  methods  following  the 
foregoing  revision  are  presented  In  Tables  3.1  and  3.2.  Except  for  the  case 
with  NCELL  -  15  and  OUTFR  ■■  0.99,  the  virtual  cell  analysis  does  very  well. 
The  characteristic  based  method  actually  shows  a  larger  disparity  from  the 
analytical  prediction  of  base  pressure.  We  note  however,  that  the  results  for 
the  physical  element  are  almost  identical  to  those  obtained  with  TDNOVA 
(Robbins,  1983);  namely;  exit  time,  10.58  ms;  velocity,  808.5  m/s;  and  base 
pressure,  54.41  HFa. 

Althou^  the  characteristic  based  method  involved  additional  computation 
near  the  boundary,  the  effect  on  run  time  was  not  large.  For  the  virtual 
element  analysis  the  CPU  times  on  a  Silicon  Graphics  Personal  Iris  workstation 
were  0.5,  2.0  and  7.3  seconds  for  15,  30  and  60  cells,  while  for  the 
characteristic  based  method  the  times  were  0.6,  2.0  and  7.3  seconds. 

Based  on  the  results  of  Tables  3.1  and  3.2  there  would  seem  to  be  no 
reason  not  to  use  the  virtual  cell  analysis.  However,  inspection  of  the 
solution  at  earlier  times  shows  that  this  method,  even  with  the  value  of  Vj, 
overwritten,  may  still  produce  significant  errors.  We  refer  to  Table  3.3. 
which  presents  results  at  0.4772  ms.  At  this  early  time  the  characteristic 
based  method  is  in  close  agreement  with  the  analytical  solution  while  the 
virtual  element  analysis  exhibits  errors  as  large  as  7.6Z,  and  these  errors  do 
not  decrease  as  NCELL  is  increased.  The  fact  that  the  final  values  are  as  good 
as  they  are  is  a  tribute  to  the  stabilizing  power  of  the  LCPFCT  algorithm. 
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TabI*  3.1.  Caosparlaon  of  LCPFCT  solution  with  virtual  call  aobaddlns  to  that  of  Lova  and  Pldduck 
(1921)  for  Lagransa  gun  problan  at  oaixzla  axlt.  WCELL  la  tha  nunbaz  of  calls.  OUTFR  la  tha 
fraction  of  tha  last  call  In  front  of  tha  projactlla  baaa,  Tha  Courant  ntanbar  la  0.5. 


Method 

HCELL 

OUTFR 

Exit  Tlioa 
(msac) 

X  Dlff 

Velocity 

(m/a) 

X  Dlff 

Base  Pressure  X  Dlff 
(MPa) 

Analytic 

- 

- 

10.56 

- 

607.7 

- 

56.19 

-  1 

LCPFCT 

15 

0.00 

10.55 

“0.28 

609.8 

0.26 

56.63 

0.66 

0.25 

10.56 

-0.19 

809.6 

0.21 

56.62 

0.62 

0.50 

10.56 

-0.19 

609.2 

0.19 

56.62 

0.42 

0.75 

10.56 

-0.19 

609.1 

0.17 

56.29 

0.16 

0.90 

10.57 

-0.10 

606.9 

0.15 

56.63 

0.66 

0.99 

10.56 

0.00 

808.3 

0.07 

55.33 

2.10 

LCPFCT 

30 

0.00 

10.56 

-0.19 

809.2 

0.19 

56.37 

0.33 

0.25 

10.57 

-0.10 

609.0 

0.16 

56.36 

0.35 

0.50 

10.57 

1 

o 

M 

O 

806.9 

0.15 

56.39 

0.37 

0.75 

10.57 

-0.10 

606.6 

0,16 

56.35 

0.30 

0.90 

10.57 

1 

o 

»-• 

o 

608.7 

0.12 

56,31 

0.22 

0.99 

10.56 

0.00 

606.5 

0.10 

54.05 

-0.26 

LCPFCT 

60 

0.00 

10,57 

-0.10 

808.9 

0.15 

56.61 

0.61 

0.25 

10.57 

-0.10 

606.6 

0.16 

56.39 

0.37 

0.50 

10.57 

-0.10 

606.7 

0.12 

56.60 

0.39 

0.75 

10.57 

o 

o 

1 

606.7 

0.12 

56.39 

0.37 

0.90 

10.56 

0.00 

606.6 

0.11 

56.36 

0.31 

0.99 

10.56 

0.00 

808.5 

0.10 

56.07 

-0.22 

Table  3.2.  Cooparlson  of  LCPFCT  solution  with  characteristic  analysis  of  boundary  cell  to  that  of 

Love  and  Pidduck  (1921)  for  Lagrange  gun  problem  at  muzzle  exit.  HCELL  is  the  number  of  calls. 

OUTFR  is  the  fraction  of  the  last  cell  in  front  of  the  projectile  base.  The  Courant  number  la  O.S. 

Method 

KCELL 

OUTFR 

Exit  Time 
(msec) 

X  Diff 

Velocity 

(m/s) 

X  Diff 

Base  Pressure  X  Diff 
(MPa) 

Analytic 

- 

- 

10.58 

- 

807.7 

- 

54,19 

- 

LCPFCT 

15 

0.00 

10.56 

0.00 

606.5 

0.10 

54.55 

0.66 

0.25 

10.57 

-0.10 

806.6 

0.11 

56.55 

0.66 

0.50 

10.56 

0.00 

806.6 

0.09 

56.69 

0.55 

0.75 

10.56 

0.00 

806.3 

0.07 

5«.42 

0.62 

0.90 

10.58 

0.00 

608.3 

0.07 

56.37 

0.33 

0.99 

10.58 

0.00 

606.3 

0.07 

56.35 

0.30 

LCPFCT 

30 

0.00 

10.57 

o 

o 

1 

606.6 

0.16 

56.51 

0.59 

0.25 

10.57 

-0.10 

606.7 

0.12 

56.50 

0.57 

0.50 

10.56 

0.00 

608.6 

0.11 

56.68 

0.56 

0.75 

10.58 

0.00 

606.5 

0.10 

56.45 

0.46 

0.90 

10.56 

0.00 

606.5 

0.10 

56.63 

0.66 

0.99 

10.56 

0.00 

606.5 

0.10 

56.61 

0.61 

LCPFCT 

60 

0.00 

10.56 

0.00 

606.7 

0.12 

56.66 

0.50 

0.25 

10.58 

0.00 

608.7 

0.12 

56.65 

0.68 

0.50 

10.58 

0.00 

608.7 

0.12 

56.66 

0.66 

0.75 

10.58 

0.00 

808.6 

0.11 

56.62 

0.62 

0.90 

10.56 

0.00 

608.6 

0.11 

54.60 

0.39 

0.99 

10.56 

0.00 

606.5 

0.10 

56.65 

0.68 
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Table  3.3.  Cooipariaon  of  base  pressures  at  0.4772  ds  ccxnputed  usins 
analytic  setbod  of  Love  and  Pidduck  (1921)  witb  LCFFCT  with  virtual 
cell  (VC)  and  LCFFCT  with  characteristic  OMthod  (CM).  The  fraction 
of  the  last  call  outside  the  boundary  is  0.99  for  both  the  VC  and  CM 
calculations.  The  Courant  nuzober  is  0.5. 


Method 

HCELL 

Base  Praasura 
(MPa) 

X  Dlff 

Analytic 

- 

554.2 

- 

LCPFCT/VC 

15 

582.1 

5.03 

30 

544.1 

-1.62 

60 

512.2 

-7.58 

LCPFCT/CM 

15 

554.3 

0.02 

30 

554.2 

0.00 

60 

554.1 

-0.02 

On  balance,  it  appears  that  the  virtual  cell  analysis  could  be  a 
satisfactory  method  under  some  circumstances,  particularly  if  we  treat  cells 
which  are  almost  completely  obstructed  as  closed.  However,  the  characteristic 
based  method  does  seem  to  offer  greater  accuracy  and  stability  without  much 
additional  cost.  Of  course,  in  a  highly  parallel  computing  environment,  the 
additional  work  at  the  boundary  will  slow  down  the  entire  calculation  at  a 
greater  proportional  rate  than  in  the  present  serial  calculation.  Yet,  if  the 
method  allows  the  same  degree  of  accuracy  to  be  achieved  with  fewer  points, 
considerable  savings  in  run  time  can  be  achieved,  especially  in  three- 
dimensional  calculations. 

We  have  selected  the  characteristic  based  analysis  for  the  following 
reasons.  Although  the  present  study  shows  the  virtual  cell  analysis  to  be 
satisfactory  for  all  but  the  most  strongly  occluded  cells  we  believe  that 
practical  limits  on  mesh  size  make  it  desirable  to  maintain  stability  for  cells 
which  are  as  much  as  99i  occluded.  Especially  when  cells  are  located  near  an 
outer  radial  boundary,  they  can  represent  a  significant  fraction  of  the  total 
volume.  Moreover,  the  present  study  has  addressed  cells  whose  geometry  is 
essentially  stationary.  In  a  two-dimensional  application  we  expect  cells  to 
grow  or  shrink  due  to  relative  motion  of  the  external  boundaries.  Not  only  is 
it  desirable  to  follow  the  history  of  the  flow  in  a  continuous  manner  as 
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certain  cells  shrink  to  nothing  or  expand  from  nothing,  we  expect  additional 
numerical  strain  to  follow  from  the  change  In  volume  as  we  approach  zero.  This 
latter  concern  may  be  understood  from  an  examination  of  Equation  3.1.4. 
Finally,  we  look  beyond  the  present  application  to  the  computational  challenge 
presented  by  the  analysis  of  increment  containers.  It  Is  expected  that  a 
characteristic  based  method  will  be  of  great  value  In  determining  the  coupled 
Internal  boundary  values  corresponding  to  the  flow  on  each  side  of  the 
container. 

3.3  Characteristic  Based  Application  of  LCPFCT  to  Two-Dimensional  Flow  with 
Non-Uniform  Boundaries 


We  begin  by  discussing  the  representation  of  gun  and  projectile  geometry. 
Figure  3.2  illustrates  the  embedding  of  tube  and  afterbody  geometry  Into  a 
structured  rectangular  mesh  as  Is  presently  done  in  the  NGEN  Code.  It  Is 
assumed  that  the  breechface  and  the  base  of  the  projectlle/sabot  are  flat  and 
so  Is  the  base  of  the  afterbody  If  one  Is  present.  The  tube  and  afterbody  are 
allowed  to  have  arbitrary  geometry  as  expressed  by  tables  of  radii  versus  axial 
distance.  It  Is  assumed  that  the  structure  of  the  geometry  Is  consistent  with 
the  resolving  power  of  the  mesh.  Also,  for  the  time  being,  we  constrain  the 
geometry  so  that  each  axial  or  radial  sweep  with  LCPFCT  Involves  a  simply 
connected  set  of  cells.  This  effectively  limits  the  afterbody  shape  to  one  in 
which  the  radius  does  not  decrease  with  distance  down  the  tube.  This  restric¬ 
tion  is  in  no  way  fundamental  and  can  be  removed  by  the  extension  of  the  coding 
to  look  for  multiple  boundary  intersections  on  axial  and  radial  lines. 

The  axial  distribution  of  cells  is  such  that  a  fixed  uniform  complement 
is  spaced  between  the  breechface  and  the  base  of  the  afterbody,  if  present,  and 
expands  in  concert  with  the  motion  of  the  projectile.  Another  fixed  uniform 
complement  is  attached  to  the  afterbody,  if  present,  and  moves  with  the 
projectile.  The  radial  distribution  is  divided  into  two  regions.  The  cells 
in  the  inner  region  are  evenly  spaced  between  the  centerline  and  the  outside 
of  the  base  of  the  afterbody.  The  cells  in  the  outer  region  are  evenly  spaced 
between  the  outside  of  the  base  of  the  afterbody  and  the  largest  radius  of  the 
chamber . 
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If  no  afterbody  is  present  the  distribution  is  axially  uniform  from  the 
breech  to  the  base  of  the  projectile  and  radially  uniform  from  the  centerline 
to  the  largest  chamber  radius.  When  an  afterbody  is  present,  motion  of  the 
projectile  will  eventually  result  in  a  large  difference  in  axial  cell  spacing 
at  the  base  of  the  afterbody.  This  could  be  avoided  by  the  Introduction  of 
coding  to  increase  the  number  of  cells  behind  the  afterbody  from  time  to  time. 

At  each  time  step  the  radial  boundaries  are  scanned  for  intersections 
with  the  array  of  cells.  Cells  are  marked  as  open,  partially  occluded  or 
closed.  Closed  cells  are  ignored  in  the  calculation,  but  they  are  assigned 
default  state  variable  values  in  order  to  facilitate  plotting.  Partially 
occluded  cells  are  assigned  values  of  volume  and  area  fractions  which  are 
transmitted  to  the  LCPFCT  metric  routines  to  modify  the  nominal  values 
corresponding  to  open  cells.  The  volume  and  area  fractions  are  computed 
analytically  on  the  assumption  that  the  boundary  element  is  linear  across  the 
cell  in  question. 

Because  the  balance  equations  are  solved  in  divergence  form,  it  is 
necessary  to  impose  a  geometric  conservation  law  (Thomas  and  Lombard,  1979). 
The  volume  fraction  is  only  computed  analytically  at  the  first  step.  Subse¬ 
quently  it  is  deduced  from  an  integral  of  the  area  fraction  times  the  normal 
cell  face  velocity.  This  eliminates  the  Introduction  of  a  fictitious  source 
term  due  to  inconsistencies  between  the  finite  volume  and  area  fractions. 

We  now  indicate  how  the  characteristic  analysis  is  applied  to  the 
multidimensional  situation  with  non-uniform  boundaries.  Cells  are  first 
integrated  in  the  axial  direction  at  all  radial  locations.  The  integration  is 
performed  from  the  first  open  or  partially  occluded  cell  at  each  radial 
position  to  the  last  such  cell,  it  being  assumed,  as  noted  above,  that  only  one 
such  string  exists.  Boundary  cells  adjacent  to  the  flat  breechface  or 
proJectile/sabot  base  are  always  open  so  that  their  treatment  is  straight¬ 
forward.  Boundary  cells  created  by  termination  by  the  tube  wall  or  the 
afterbody  may  be  partially  occluded.  They  are  treated  after  the  fashion  of  the 
virtual  cell  analysis  described  in  the  previous  section.  Thus,  on  the  axial 
sweep,  we  only  account  for  area  and  volume  fraction. 
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When  the  axial  sweep  is  cooplete  at  all  locations  we  perform  the  radial 
sweep  at  each  axial  location.  It  is  at  this  point  that  the  characteristic 
analysis  is  used  to  stabilize  the  solution.  We  use  the  external  boundary 
position  to  define  the  boundaries  of  the  first  and  last  cells  to  be  Integrated 
and  we  Introduce  the  normal  to  the  boundary  in  setting  the  boundary  values . 
The  axial  velocity  component  is  extrapolated  to  the  wall  and  the  value  of 
radial  velocity  at  the  wall  is  determined  from  the  condition  of  attachment. 
This  value  is  then  used  together  with  the  interior  data  to  compute  a  character¬ 
istic  value  of  bovmdary  pressure  which  is  then  compared  with  the  value  implied 
by  the  L^PFCT  integration  as  in  the  one -dimensional  case  discussed  previously. 
Again,  mass  transfer  between  the  boundary  cell  and  its  neighbor  is  introduced 
to  bring  the  two  values  into  agreement.  The  influence  of  the  cross  flow  terms 
is  taken  to  be  embedded  in  the  LCPFCT  axial  sweep  so  that  the  characteristic 
analysis  is  effectively  one -dimensional  as  far  as  spatial  differencing  is 
concerned,  thereby  minimizing  the  computational  burden. 


Appendix  A  discusses  the  characteristic  analysis  of  a  two-dimensional 
single  phase  flow  with  a  non-homogeneous  source  term.  The  results  contained 
therein  form  the  basis  for  the  present  application  since  we  formally  include 
with  the  non-homogenous  terms  the  derivative  of  porosity.  Equation  A. 2. 27  is 
the  relevant  characteristic  form,  and  we  restate  it  here  as 
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Now  includes  the  non-homogenous  terms  and  the  crossflow  (axial)  derivatives. 
We  have  v,,  as  the  radial  mesh  velocity  and  is  defined  by 
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where  the  ,  Including  which  is  not  coupled,  are  defined  as  follows 
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We  write  the  as  ^i,  +  f it  t  i  -  1  -  to  conform  with  the  timesplitting 
convention.  Then  we  have  (Gough,  1983),  with  the  formal  representation  of  the 
porosity  derivative  as  non-homogeneous 
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3.3.13 


^41  =  0  .  3.3.14 

Now  the  2-conrponents  of  can  be  resolved  as  the  result  of  the  integration  of 
the  axial  equations  of  motion  including  the  relevant  source  terms.  This  is 
exactly  what  occurs  on  the  axial  sweep  by  LCPFCT  except  that  the  influence  of 
changing  cell  volume  is  also  captured.  Then,  factoring  out  this  term  and  using 
a  superscript  ~  to  denote  the  result  of  the  axial  sweep,  Equation  3.3.1  may  be 
discretized  at  the  boundary  in  the  form 


where  involves  only  the  . 


The  spatial  derivatives  are  understood  to  be  evaluated  using  old  data. 
Ve  understand  DA/Dt  to  be  the  axial  convective  derivative  of  the  cell  cross- 
sectional  area.  This  is  the  only  additional  term  which  needs  to  be 
differenced.  Equation  3.3.15  therefore  makes  only  minimal  computational 
demands . 

3.4  The  Larpe  Particle  Integrator  (LPI) 

Ue  now  discuss  the  LPI  which  is  intended  to  be  applied  to  those  discrete 
phases  which  cannot  be  assumed  to  be  in  mechanical  equilibrium  with  the 
continuous  phase.  In  all  cases  of  practical  interest  this  will  include  all 
solid  propellant  increments. 

In  previous  work  (Gough,  1983)  we  have  explained  the  importance  of  the 
presence  of  ullage  in  solid-propellant  charges  and  of  the  properties  of  the 
increment  containers  which  are  generally  present  in  artillery  charges.  The 
initial  presence  and  subsequent  persistence  of  ullage  and  the  permeability  and 
fracture  properties  of  the  containers  can  exert  considerable  influence  over  the 
path  of  flamespreading  and  the  subsequent  ballistic  behavior  of  the  charge. 
Small  differences  in  these  attributes  can  make  the  difference  between  a  safe 
charge  and  one  which  can  destroy  the  gun  (Hay  and  Horst,  1979).  Therefore,  in 
previous  work  we  have  taken  a  modeling  approach  based  on  the  representation  of 
charge  Increment  boundaries  as  explicit  discontinuities  across  which  not  only 
the  porosity  but  also  the  gas-phase  properties  may  Jump  discontlnuously .  The 
increment  boundaries  were  tracked  explicitly  as  part  of  the  numerical  solution 
process.  The  internal  boundary  conditions  linking  the  state  variables  on  each 
side  consisted  of  finite  balances  of  mass,  momentum  and  energy.  By  viewing  the 
container  as  a  surface  attribute  of  the  increment  ve  were  able  to  embed 
properties  of  reactivity  and  permeability  into  the  finite  balance  equations  or 
Jump  conditions.  Thus  our  earlier  approach  not  only  assured  precise  tracking 
of  ullage,  with  a  complete  elimination  of  the  possibility  of  contamination 
through  numerical  diffusion  of  the  solid  propellant  properties  outside  the 
increment  boundary,  but  it  also  admitted  a  representation  of  the  container 
characteristics  on  the  assumption  that  the  motion  of  the  container  was  tied  to 
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that  of  the  propellant.  Although  this  approach  was  successful,  It  Involved  a 
great  deal  of  complex  code  and  grid  management  and  a  simpler  approach  was 
sought  for  NGEN, 

There  Is  no  escaping  the  complexities  associated  with  the  presence  of 
containers.  They  must  be  treated  as  Internal  boundaries.  The  approach  taken 
here  in  respect  to  the  modeling  of  the  continuous  phase  is  thought  to 
anticipate  this  eventual  goal  in  an  appropriate  manner.  At  the  same  time,  the 
choice  of  continuum  flow  solver  has  been  made  with  a  view  to  the  simplification 
of  the  representation  of  the  increment  boundaries.  These  are  still  tracked 
explicitly  via  the  adoption  of  a  Lagrangian  representation  of  each  increment. 
But  the  boundaries  are  not  viewed  as  discontinuities.  Instead,  the  porosity 
is  made  to  vary  continuously  over  a  short  distance.  Internal  jvunp  conditions 
are  not  applied.  The  continuum  flow  solver  is  required  to  be  robust  enough  to 
integrate  the  flow  In  a  stable  manner  in  the  presence  of  strong  porosity 
gradients.  We  have  previously  shown  that  LCPFCT  meets  this  requirement  (Gough, 
1992). 


Consider  a  single  charge  increment.  In  the  present  version  of  the  code 
Its  Initial  distribution  is  defined  by  a  rear  and  forward  delimiter  and  by  an 
inner  and  an  outer  delimiter.  The  increment  Is  assumed  to  occupy  uniformly  the 
cylindrical  domain  defined  by  these  delimiters  minus  any  intrusions  by  the 
chamber  and  afterbody.  The  increment  is  then  represented  by  a  structured 
array  of  LPI -particles .  These  are  arranged  in  an  axially  uniform  manner  from 
the  rear  to  the  forward  delimiter  and  radially  from  the  inner  to  the  outer 
delimiter.  LPI-particles  which  would  lie  outside  the  tube  wall  or  inside  the 
afterbody  according  to  this  prescription  are  then  pushed  radially  to  the 
appropriate  radial  boundaries.  The  redundant  particles  are  then  assigned  zero 
weight  and  play  no  part  in  the  subsequent  solution  process. 


Each  particle  is  assigned  a  weight  as  follows.  The  grain  nuzaber  density 
may  be  computed  from 


n 


o 


Mj 

Wi 
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where  Is  the  Increment  mass  and  the  Increment  volume  and  Pp,  Vp,  are 
respectively  the  density  of  the  propellant  and  the  Initial  volume  of  an 
Individual  grain.  Then  the  number  weighting  of  a  particle  Is  taken  to  be  n^ 
times  the  volume  of  a  cylinder  defined  by  the  midpoints  between  the  particle 
and  Its  four  nel^bors.  For  boundary  particles  the  relevant  midpoint  value  Is 
replaced  by  the  appropriate  axial  or  radial  coordinate  of  the  particle.  It  Is 
easy  to  see  that  this  process  automatically  assigns  non-zero  weight  to  at  most 
one  particle  which  was  pushed  Into  contact  with  a  radial  boundary. 


In  addition  to  number  weighting,  each  LFI-partlcle  is  assigned  the 
following  attributes;  axial  and  radial  position  and  velocity,  surface 
regression,  surface  temperature  and  cubic  profile  thermal  parameter  (Equation 
2. 3. 6. 2),  and  values  of  porosity  and  Intergranular  stress. 

At  each  time  step  the  equations  of  motion  are  Integrated  using  a  simple 
first  order  time  differencing  scheme.  The  positions  are  updated  first  for  all 
LPI -particles .  Then,  In  a  second  sweep  the  velocities  are  updated  and 
subjected  to  the  external  boundary  condition.  Depending  on  the  surface 
temperature  either  the  thermal  parameter  or  the  surface  regression  is 
Integrated.  In  the  same  sweep  all  the  Interphase  transfer  properties  are 
computed  for  each  particle  and  mapped  onto  the  grid  for  the  continuum  flow 
solver  together  with  the  porosity  distributions. 

Ue  now  amplify  on  the  previous  paragraph.  First,  with  regard  to  the 
equation  of  motion,  we  see  from  an  inspection  of  Equation  2.2.14  that  we 
require  gradients  of  pressure  and  of  intergranular  stress.  These  are  computed 
quite  differently.  The  pressure  gradient  Is  computed  on  the  continuum  mesh  by 
means  of  central  differencing  to  define  cell  centered  values.  Then  for  each 
LPI -particle  a  local  value  of  pressure  gradient,  together  with  all  other 
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necessary  continuous  phase  state  variables,  is  computed  using  linear  interpola¬ 
tion  In  both  the  radial  and  the  axial  directions.  The  gradient  of  intergranu¬ 
lar  stress  Is  computed  using  differences  between  the  LPI -particles  which 
neighbor  the  particle  In  question.  One-sided  differences  -are  used  at 
boundaries . 


Ve  now  explain  the  mapping  of  LPI  data  back  to  the  continuum  mesh.  At 
each  time  step  each  particle  is  assigned  a  rectangular  domain  of  influence  as 
follows.  Let  Zi  and  be  the  left  and  rl^t  axial  delimiters  and  r^  and  r,,  be 
the  Inner  and  outer  delimiters.  Then  is  set  equal  to  the  axial  position  of 
the  particle  to  the  rear,  provided  that  such  a  particle  exists.  If  we  are  at 
the  boundary  or  if  the  neighbor  in  question  has  zero  weight,  the  is  set 
equal  to  a  nominal  value  defined  by  the  user.  If  this  value  Is  less  than  the 
position  of  the  rear  boundary  it  is  replaced  by  the  position  of  the  boundary. 
Then  z^,  r^  and  r^  are  all  defined  in  an  analogous  manner.  Let  the  coordinates 
of  the  particle  be  (z^j,  r^j )  and  let 

Vij  “  -J(Zr  -  2i,)[(rij  *  -  (r^j  +  3.4.2 

O 

be  the  volume  associated  with  the  particle.  This  is  identical  with  the  pro¬ 
cedure  used  to  define  the  number  weighting,  except  at  the  boundaries  where  the 
definition  of  the  delimiter  admits  a  deliberate  extension  into  the  ullage  for 
the  purpose  of  smoothing. 


If  respectively  the  grain  volume  and  the  number  weight¬ 

ing,  we  can  define  the  particle  volume  fraction 
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and  similarly  for  the  other  interphase  data  such  as  mass  and  enthalpy  transfer 
and  drag.  The  value  of  so  defined  is  used  to  update  the  solid  propellant 
intergranular  stress  by  means  of  a  simple  first  order  time  difference. 
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To  assign  the  LPI*particle  volume  to  the  continuum  mesh  ve  simply  scan 
the  continuum  mesh  to  find  those  cells  whose  centers  are  in  the  domain  of 
influence.  Let  (z^.r^,)  be  the  coordinates  of  such  a  cell  center.  Define 


(2*  "  2c)/(2r  -  Zij) 

4 

(Ze  -  2L)/(Zij  -  Zl)  ifZe<Zij 
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and  define  analogously  for  the  radial  direction.  Then  the  continuum  cell 
in  question  is  assigned  a  volume  fraction  contribution  and  similarly  for 
the  other  interphase  properties.  These  values  are  summed  over  all  LPI- 
particles.  It  is  not  hard  to  see  that  this  is  tantamount  to  a  bidirectional 
linear  interpolation  process. 


Some  difficulties  can  arise  with  this  approach  when  the  LPI -particle  is 
close  to  an  external  boundary  whose  shape  does  not  conform  well  with  the 
rectangular  domain  of  influence.  Local  pockets  of  elevated  porosity  can  be 
produced.  This  problem  is  expected  to  be  remedied  through  minor  finessing  of 
the  logic  to  determine  the  boundaries  of  the  domain  of  influence  in  such  cases. 


A  more  serious  potential  problem  relates  to  the  behavior  of  the  aggregate 
as  a  whole.  The  presence  of  the  pressure  gradient  in  the  solid  propellant 
equation  of  motion  can  create  a  Helmholtz  type  of  instability.  Suppose  that 
we  have  a  local  minimum  of  porosity  in  a  macroscopically  one -dimensional  flow. 
Porosity  affects  the  gas-phase  properties  in  the  same  manner  as  cross-sectional 
area  in  a  duct.  Thus,  for  subsonic  flow,  the  porosity  minimum  will  be 
accompanied  by  a  pressure  minimum.  The  pressure  gradient  will  therefore  act 
to  drive  particles  towards  the  location  of  the  minimum  of  porosity,  deepening 
it.  The  process  may  possibly  be  stabilized  to  some  extent  by  interphase  drag 
and  certainly  by  Intergranular  stresses  as  the  porosity  becomes  small  enough. 
Nevertheless,  strong  variations  in  porosity  can  and  will  arise  and  are  Inherent 
in  the  formulation  of  the  equations.  It  is  not  presently  known  how  best  to 
address  this  problem. 
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To  assign  the  LPI •particle  volume  to  the  continuum  mesh  we  simply  scan 
the  continuum  mesh  to  find  those  cells  whose  centers  are  In  the  domain  of 
Influence.  Let  (Ze,re)  be  the  coordinates  of  such  a  cell  center.  Define 


w. 
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(Ze  -  Zi,)/(2ij  -  2l) 


if  Ze  ^  Zij 
if  Ze  <  Zij 
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and  define  analogously  for  the  radial  direction.  Then  the  continuum  cell 
In  question  Is  assigned  a  volume  fraction  contribution  w^w^a^^j  and  similarly  for 
the  other  Interphase  properties.  These  values  are  summed  over  all  LPI- 
partlcles.  It  is  not  hard  to  see  that  this  is  tantamount  to  a  bidirectional 
linear  Interpolation  process. 


Some  difficulties  can  arise  with  this  approach  when  the  LFl-particle  Is 
close  to  an  external  boundary  whose  shape  does  not  conform  well  with  the 
rectangular  domain  of  Influence.  Local  pockets  of  elevated  porosity  can  be 
produced.  This  problem  Is  expected  to  be  remedied  through  minor  finessing  of 
the  logic  to  determine  the  boundaries  of  the  domain  of  Influence  In  such  cases. 


A  more  serious  potential  problem  relates  to  the  behavior  of  the  aggregate 
as  a  whole.  The  presence  of  the  pressure  gradient  In  the  solid  propellant 
equation  of  motion  can  create  a  Helmholtz  type  of  Instability.  Suppose  that 
we  have  a  local  minimum  of  porosity  In  a  macroscoplcally  one -dimensional  flow. 
Porosity  affects  the  gas -phase  properties  In  the  same  manner  as  cross-sectional 
area  In  a  duct.  Thus,  for  subsonic  flow,  the  porosity  minimum  will  be 
accompanied  by  a  pressure  minimum.  The  pressure  gradient  will  therefore  act 
to  drive  particles  towards  the  location  of  the  minimum  of  porosity,  deepening 
It.  The  process  may  possibly  be  stabilized  to  some  extent  by  Interphase  drag 
and  certainly  by  Intergranular  stresses  as  the  porosity  becomes  small  enough. 
Nevertheless,  strong  variations  In  porosity  can  and  will  arise  and  are  Inherent 
in  the  formulation  of  the  equations.  It  Is  not  presently  known  how  best  to 
address  this  problem. 
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4.0  NUMERICAL  SOLUTIONS 

The  operability  of  the  NGEN  Code  and  Its  general  applicability  to  solid 
propellant  charges  are  illustrated  here  by  reference  to  two  different  data 
bases.  The  first  of  these  represents  a  ISSmn  howitzer  charge  and  the  second 
represents  a  120inm  tank  gun  charge.  In  both  cases  we  represent  the  actual 
geometry  of  the  tube  to  an  extent  consistent  with  the  capabilities  of  the  XKTC 
Code.  Thus  it  is  assumed  that  the  breech  face  and  the  base  of  the  projectile 
(or  the  base  of  the  sabot)  are  flat.  The  tube  wall  geometry  is  captured  as  is 
the  geometry  of  the  afterbody  in  the  case  of  the  tank  gun  simulation. 

ISSmm  Howitzer  Simulation 

The  complete  NGEN  input  data  file  for  this  problem  is  presented  in  Table 
4.1.  The  data  base  was  developed  from  the  one ‘dimensional  XKTC  (Gough,  1986) 
data  base  for  the  M203  propelling  charge.  For  simplicity  flamespreading  is  not 
considered.  The  charge  is  taken  to  be  ignited  at  the  initial  instant.  The 
discharge  from  the  igniter  is  ignored  and  the  initial  condition  corresponds  to 
ambient  temperature  and  pressure.  Bore  resistance  is  likewise  ignored.  The 
propellant  consists  of  seven-perforation  grains  with  the  correct  geometry  and 
thermodynamic  properties.  However,  the  bum  rate  is  adjusted  to  give 
approximately  the  computed  maximum  pressure  for  the  original  XKTC  data  base. 

The  chamber  radius  is  taken  to  taper  linearly  with  axial  distance  from 
a  value  of  9.17  cm  at  the  breechface  to  a  value  of  7.82  cm  at  the  entrance  of 
the  tube.  In  the  NGEN  simulation  the  propellant  is  initially  configured  with 
ullage  in  front  of  it  and  around  it.  The  forward  boundary  delimiter  is  75  cm 
which  is  less  than  the  position  of  the  projectile  base  at  82.14  cm.  Similarly, 
the  outer  radial  delimiter  is  7.6  cm  which  is  smaller  than  the  tube  radius  at 
all  locations.  Of  course,  the  representation  of  the  charge  properties  will 
involve  some  smearing  of  the  porosity  outside  the  delimiters  as  discussed  in 
Chapter  3.0. 
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Table  4.1.  Input  Data  for  NGEN  Simulation  of  ISSmm  Howitzer 
Propelling  Charge 


CONTROL  PARAMETERS 


NPRINT(0-N0  PRINT, 1-PRINT)  1 
NSUMRY(0-N0  SUMMARY  TABLES , 1-YES )  1 
NPL0T(0-N0  ISOMETRIC  CARPET  PLOTS , 1-PLOT)  1 
NDSKW(0-N0  DISC  SAVE, 1-DISC  SAVE)  0 
NDSKR(0-N0  DISC  START, >0-DISC  START  AT  STEP 

NDSKR)  0 


ISOMETRICALLY  PLOTTED  QUANTITIES  ( 1-YES, O-NO) 

MESH  0  POROSITY  1  GRANULAR  STRESS  0  PRESSURE  1 
DENSITY  0  GAS  AXIAL  VELOCITY  1  SOLID  AXIAL  VELOCITY  1 
GAS  RADIAL  VELOCITY  1  SOLID  RADIAL  VELOCITY  1  GAS  TEMPERATURE  0 
PARTICLE  SURFACE  TEMPERATURE  0 


TERMINATION  AND  LOGOUT  PARAMETERS 


MAX.  STEPS  10000 

MAX.  TRAVEL(CM)  520.700 

STEPS  TO  LOGOUT  5000 

INTERVAL  TO  LOGOUT(MSEC)  0.500 

DEBUG  PRINT ( O-NO ; 1-YES )  0 


INTEGRATION  DATA 


COURANT  NUMBER(-)  0.500 

ANTIDIFFUSIVE  FLUX  MULTIPLIER  0.999 

TOTAL  NUMBER  OF  RADIAL  CELLS(-)  8 

MINIMUM  NUMBER  OF  AXIAL  CELLS  BEHIND  AFTERBODY  30 

MAXIMUM  NUMBER  OF  AXIAL  CELLS  BEHIND  AFTERBODY  30 

NUMBER  OF  AXIAL  CELLS  ALONG  AFTERBODY  0 

NUMBER  OF  RADIAL  CELLS  AT  BASE  OF  AFTERBODY  0 


INITIAL  DATA 


PRESSURE(MPA)  0.101 
TEMPERATURE (K)  294.0 
GAMMA(-)  1.2430 
MOL . WGT ( GM/GMOL)  2  3 . 4400 
C0V0LUME(CM**3/GM)  1.0300 
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CHARGE  REPRESENTATION  PARAMETERS 


NUMBER  OF  PROPELLANTS  IN  CHARGE  1 
NUMBER  OF  INCREMENTS  IN  CHARGE  1 
PROPELLANT  MODEL  (1-LP,2-2D)  2 


PROPERTIES  OF  PROPELLANT  TYPE  1 


SOLID  PHASE  CONSTITUTIVE  DATA 


SETTLING  POROSITY  OF  GRANULAR  BED(-) 
SPEED  OF  COMPRESSION  WAVE(CM/SEC) 
SPEED  OF  EXPANSION  WAVE(CM/SEC) 
DENSITY  OF  SOLID  PHASE(GM/CC) 
THERMAL  CONDUCTIVITY(J/CM-SEC-DEG.K) 
THERMAL  DIFFUSIVITY(CM**2/SEC) 


0.40000 

15240.0 

127000.0 

1.5830 

0.160100E-02 

0.645160E-03 


GAS  PHASE  CONSTITUTIVE  DATA  AND  PROPELLANT  INITIAL 


RATIO  OF  SPECIFIC  HEATS(-)  1.24300 
MOLECULAR  WEIGHT(GM/GM-MOL)  23.440 
COVOLUME (CC/GM)  1.030 
INITIAL  TEMPERATURE  (DEG.K)  445.0 


SOLID  PHASE  COMBUSTION  CHARACTERISTICS 

2 

444 . 4 
4384. 


NUMBER  OF  BURN  RATE  DATA 
IGNITION  TEMPERATURE (DEG.K) 
CHEMICAL  ENERGY(J/GM) 


MAX.  PRESSURE 
(MPA) 


BURN  RATE  ADDITIVE  CONSTANT  PRE- EXPONENT 
(CM/SEC)  (CM/SEC-MPA**BN) 


68.95 

689.50 


0.00000 

0.00000 


0.24810 

0.24810 


GRAIN  GEOMETRY 


FORM(0-CYLINDER, 1-SPHERE, 2-SOLID  STICK,  0 
3-PERFORATED  OR  SLOTTED  STICK) 

EXTERNAL  DIAMETER(CM)  1.060 
LENGTH (CM)  2.408 
DIAMETER  OF  PERFORATIONS (CM)  0.086 
NUMBER  OF  PERFORATIONS(-)  7. 


TEMPERATURE 


EXPONENT 

(-) 

0.78640 

0.78640 
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PROPERTIES  OF  INCREMENT  NUMBER 


1 


MAIN  CHARGE  PROPELLANT  TYPE 
NUMBER  OF  AXIAL  CELLS 
NUMBER  OF  RADIAL  CELLS 
MASS  OF  MAIN  CHARGE (G) 

REAR  BOUNDARY  DELIMITER  (CM) 
FORWARD  BOUNDARY  DELIMITER  (CM) 
INNER  BOUNDARY  DELIMITER  (CM) 
OUTER  BOUNDARY  DELIMITER  (CM) 


1 

30 

8 

11861.4000 

0.000 

75.000 

0.000 

7.600 


PROPERTIES  OF  PROJECTILE 


INITIAL  POSITION  OF  PROJECTILE/SABOT  BASE (CM) 
PROJECTILE  MASS(GM) 


82.140 

43181.800 


NUMBER  OF  ENTRIES  IN  BORE  RESISTANCE  TABLE  2 

RESISTANCE  LAW  NUMBER  0 

N.B.  IF  <1  OR  >3,  VALUE  WILL  DEFAULT  TO  2  INTERNALLY 


BORE  RESISTANCE  DATA 

PROJECTILE  TRAVEL(CM)  RESISTIVE  PRESSURE(MPA) 

0.000  0.000 

600.000  0.000 

NUMBER  OF  TUBE  GEOMETRY  DATA  3 


TUBE  GEOMETRY 


AXIAL  POS(CM)  RADIUS (CM) 


O.OOOE+00 

82.1 

602. 


9.17 

7.82 

7.82 
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The  continuum  mesh  is  given  8  radial  and  30  axial  cells.  Due  to  the 
taper  of  the  chamber,  not  all  cells  are  open.  The  charge  is  represented  by  8 
radial  and  30  axial  particles.  The  maximum  Courant  number  is  0.5. 

The  distribution  of  porosity  at  four  times  are  shown  in  Figure  4.1.  The 
distribution  at  t  -  0  ms  also  illustrates  the  continuum  grid,  it  being 
understood  that  the  plotted  lines  correspond  to  cell  centers  except  on  the 
external  boundaries.  The  initial  ullage  is  apparent  from  the  porosity 
distribution.  Me  note  that  the  smearing  at  the  boundaries  results  in  a  zone 
of  sharply  varying  but  not  discontinuous  properties.  By  4.0  ms  the  charge  has 
expanded  to  fill  the  chamber.  By  8.0  ms,  which  is  close  to  the  time  of  maximum 
pressure,  the  porosity  is  essentially  uniform  except  for  the  small  region  of 
\illage  behind  the  projectile  base  which  is  an  expected  feature  of  the  solution 
since  the  propellant  grains  are  expected  to  lag  to  a  certain  extent.  By  12.0 
ms  burnout  of  the  charge  is  occurring  with  the  process  beginning  at  the  breech 
where  the  elevated  pressure  causes  more  rapid  combustion. 

Figure  4.2  presents  the  distributions  of  pressure  at  four  similar  times. 
It  should  be  noted  that  the  viewpoint  is  from  the  tube  wall.  It  is  interesting 
to  note  how  rapidly  the  pressure  equilibrates  in  the  radial  direction.  The 
distributions  are  essentially  one •dimensional  at  all  times.  Including  the 
earliest  at  0.5  ms.  This  is  true  even  though  the  gas -phase  flow  field  at  0.5 
ms,  shown  in  Figure  4.3,  is  far  from  one -dimensional.  Strong  radial  convection 
is  associated  with  the  pressurization  of  the  radial  ullage  by  the  products  of 
combustion.  This  is  of  course  strongest  near  the  breech  where  the  radial 
ullage  is  most  pronounced.  At  the  base  of  the  projectile  the  suction  caused 
by  motion  of  the  projectile  dominates  the  flow  which  is  essentially  axial. 

As  a  check  on  the  overall  accuracy  of  the  coding,  a  comparison  was  made 
of  the  NGEN  ballistic  predictions  with  those  of  XKTC.  The  data  bases  were  made 
completely  consistent  except  for  the  representation  of  the  radial  ullage  in 
NGEN,  a  feature  absent  from  XICTC.  The  agreement  is  shown  in  Figure  4.4  and  is 
thought  to  be  satisfactory. 
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120mm  Tank  Gun  Simulations 


The  complete  NGEN  input  data  file  for  this  problem  is  presented  in  Table 
4.2.  •  The  geometrical  data  were  developed  from  an  XKTC  data  base  for  a  tank  gun 
round  with  a  strongly  intruding  afterbody.  Apart  from  the  projectile  mass,  the 
remaining  data  are  purely  nominal.  The  charge  is  represented  as  completely 
filling  the  available  initial  volume.  It  is  taken  to  be  seven  perforation 
granular  with  burn  rate  adjusted  to  give  a  moderately  high  pressure  of  nearly 
380  MPa.  Flamespreading  is  neglected,  the  charge  being  initially  ignited. 
Likewise,  we  ignore  bore  resistance  and  the  possible  presence  of  a  combustible 
cartridge  case.  The  initial  pressure  is  taken  to  be  6.895  MPa  and  the 
temperature  is  3009K. 

The  purpose  of  the  simulation  is  to  test  the  ability  of  NGEN  to  represent 
the  strongly  non-uniform  boundary  geometry  t3rpical  of  tank  gun  ammunition.  The 
continuvun  mesh  is  assigned  10  radial  cells,  with  3  assigned  to  the  base  of  the 
afterbody.  A  total  of  30  axial  cells  are  assigned  with  15  behind  the  afterbody 
and  15  along  the  afterbody.  Clearly,  due  to  the  non-uniform  geometry  many  of 
the  cells  are  closed.  The  propellant  is  represented  by  8  radial  and  30  axial 
particles.  All  the  spatial  delimiters  are  set  equal  to  zero.  As  a  default 
procedure,  NGEN  takes  the  propellant  to  occupy  the  entire  initial  volume  of  the 
chamber,  taking  into  account  the  intrusion  of  the  afterbody.  As  discussed  in 
Chapter  3.0,  those  particles  which  would  be  covered  by  the  external  boundaries 
are  assigned  zero  weight  and  are  not  updated  in  the  solution  algorithm.  The 
maximum  Courant  number  is  0.5. 

Figure  4.5  illustrates  the  porosity  distributions  at  four  times.  The 
distribution  at  t  —  0.0  ms  also  illustrates  the  continuum  grid.  It  should  be 
noted  that  the  base  of  the  afterbody  is  actually  flat.  The  shape  of  the  base 
shown  in  the  figure  is  an  artifact  of  the  plot  interface  routine.  The  figure 
illustrates  clearly  the  non-uniform  nature  of  the  boundary  geometry  and  its 
representation  by  a  grid  which  is  everywhere  rectangular  except  adjacent  to  the 
radial  boundaries.  The  distributions  show  the  gradual  separation  of  the 
propellant  from  the  base  of  the  sabot  and  the  base  of  the  afterbody.  We  also 
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Table  4.2.  Input  Data  for  NGEN  Siinulation  of  120nuii  Tank  Gun  Charge 


CONTROL  PARAMETERS 


NPRINT(0-N0  PRINT, 1-PRINT)  1 
NSUMRY(0-N0  SUMMARY  TABLES , 1-YES )  1 
NPLOT(0-N0  ISOMETRIC  CARPET  PLOTS , 1-PLOT)  1 
NDSKW(0-N0  DISC  SAVE, 1-DISC  SAVE)  0 
NDSKR(0-N0  DISC  START, >0-DISC  START  AT  STEP 

NDSKR)  0 


ISOMETRICALLY  PLOTTED  QUANTITIES  ( 1-YES, O-NO) 

MESH  0  POROSITY  1  GRANULAR  STRESS  0  PRESSURE  1 
DENSITY  0  GAS  AXIAL  VELOCITY  1  SOLID  AXIAL  VELOCITY  0 
GAS  RADIAL  VELOCITY  1  SOLID  RADIAL  VELOCITY  0  GAS  TEMPERATURE  0 
PARTICLE  SURFACE  TEMPERATURE  0 


TERMINATION  AND  LOGOUT  PARAMETERS 


MAX.  STEPS  10000 

MAX.  TRAVEL(CM)  475.000 

STEPS  TO  LOGOUT  10000 

INTERVAL  TO  LOGOUT (MSEC)  0.500 

DEBUG  PRINT (O-NO; 1-YES)  0 


INTEGRATION  DATA 


COURANT  NUMBER(-)  0.500 

ANTIDIFFUSIVE  FLUX  MULTIPLIER  0.999 

TOTAL  NUMBER  OF  RADIAL  CELLS(-)  10 

MINIMUM  NUMBER  OF  AXIAL  CELLS  BEHIND  AFTEEIBODY  15 

MAXIMUM  NUMBER  OF  AXIAL  CELLS  BEHIND  AFTERBODY  15 

NUMBER  OF  AXIAL  CELLS  ALONG  AFTERBODY  15 

NUMBER  OF  RADIAL  CELLS  AT  BASE  OF  AFTERBODY  3 


INITIAL  DATA 


PRESSURE(MPA)  6.895 
TEMPERATURE (K)  3009.0 
GAMMA(-)  1.2430 
MOL.WGT(GM/GMOL)  23.3600 
COVOLUME ( CM** 3/GM)  1.0300 
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CHARGE  REPRESENTATION  PARAMETERS 


NUMBER  OF  PROPELLANTS  IN  CHARGE  1 
NUMBER  OF  INCREMENTS  IN  CHARGE  1 
PROPELLANT  MODEL  (1-LP,2-2D)  2 


PROPERTIES  OF  PROPELLANT  TYPE  1 
SOLID  PHASE  CONSTITUTIVE  DATA 


SETTLING  POROSITY  OF  GRANULAR  BED(-) 
SPEED  OF  COMPRESSION  WAVE(CM/SEC) 
SPEED  OF  EXPANSION  WAVE(CM/SEC) 
DENSITY  OF  SOLID  PHASE(GM/CC) 

THERMAL  CONDUCTIVITY(J/CM-SEC-DEG .K) 
THERMAL  DIFFUSIVITY(CM**2/SEC) 


0.40000 

15240.0 

127000.0 

1.5830 

0.160100E-02 

0.645160E-03 


GAS  PHASE  CONSTITUTIVE  DATA  AND  PROPELLANT  INITIAL  TEMPERATURE 


RATIO  OF  SPECIFIC  HEATS(-)  1.24300 
MOLECULAR  WEIGHT (GM/GM-MOL)  23.260 
COVOLUME (CC/GM)  1.030 
INITIAL  TEMPERATURE  (DEG.K)  445.0 


SOLID  PHASE  COMBUSTION  CHARACTERISTICS 

2 

444.4 
4384. 


NUMBER  OF  BURN  RATE  DATA 
IGNITION  TEMPERATURE (DEG.K) 
CHEMICAL  ENERGY (J/GM) 


MAX.  PRESSURE  BURN  RATE  ADDITIVE  CONSTANT  PRE- EXPONENT  EXPONENT 
(MPA)  (CM/SEC)  ( CM/SEC -MPA**BN)  (-) 

68.95  0.00000  0.50000  0.70000 

689.50  0.00000  0.50000  0.70000 


GRAIN  GEOMETRY 


F0RM(0-CYLINDER ,  1~SPHERE ,  2-SOLID  STICK ,  0 
3-PERFORATED  OR  SLOTTED  STICK) 

EXTERNAL  DIAMETER(CM)  1.060 
LENGTH (CM)  2.408 
DIAMETER  OF  PERFORATIONS (CM)  0.086 
NUMBER  OF  PERFORATIONS (- )  7. 
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PROPERTIES  OF  INCREMENT  NUMBER 


1 


MAIN  CHARGE  PROPELLANT  TYPE  1 

NUMBER  OF  AXIAL  CELLS  30 

NUMBER  OF  RADIAL  CELLS  8 

MASS  OF  MAIN  CHARGE(G)  7800.0000 

REAR  BOUNDARY  DELIMITER  (CM)  0.000 

FORWARD  BOUNDARY  DELIMITER  (CM)  0.000 

INNER  BOUNDARY  DELIMITER  (CM)  0.000 

OUTER  BOUNDARY  DELIMITER  (CM)  0.000 

PROPERTIES  OF  PROJECTILE 

INITIAL  POSITION  OF  PROJECTILE/SABOT  BASE(CM)  55.900 
PROJECTILE  MASS(GM)  8936.000 

NUMBER  OF  ENTRIES  IN  BORE  RESISTANCE  TABLE  2 

RESISTANCE  LAW  NUMBER  0 


N.B.  IF  <1  OR  >3,  VALUE  WILL  DEFAULT  TO  2  INTERNALLY 


BORE  RESISTANCE  DATA 

PROJECTILE  TRAVEL(CM)  RESISTIVE  PRESSURE(MPA) 

0.000 

600.000 

NUMBER  OF  TUBE  GEOMETRY  DATA  6 


TUBE  GEOMETRY 

AXIAL  POS(CM)  RADIUS (CM) 


i.OOOE+00 

5.72 

7.62 

7.85 

48.3 

7.85 

55.9 

6.05 

59.4 

6.00 

529. 

6.00 

NUMBER  OF  AFTEEIBODY  GEOMETRY  DATA  2 


0.000 

0.000 


AFTERBODY  GEOMETRY 


AXIAL  POS(CM)  RADIUS (CM) 


0 . OOOE+00  1 . 90 

36.8  2.67 
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note,  at  4.5  ms,  the  separation  from  the  re-entrant  comer  at  the  mouth  of  the 
tube,  as  grains  are  swept  around  the  chambrage. 

Additional  details  of  the  simulation  are  shown  in  Figures  4.6  -  4.8, 
which  illustrate  the  gas-phase  flow  field  at  three  times,  and  by  Figures 
4.9  -  4.11,  which  illustrate  the  pressure  contours  at  the  same  times.  It  is 
interesting  to  note  the  suction  created  by  the  base  of  the  afterbody  which 
strongly  perturbs  the  flow  field  away  from  the  axial  structure  implicit  in  one- 
dimensional  codes.  We  note  that  since  the  grid  behind  the  afterbody  is 
expanding  in  concert  with  the  motion  of  the  projectile,  while  that  along  the 
afterbody  retains  its  initial  axial  spacing,  there  is  a  strong  discontinuity 
in  cell  size  near  the  base  of  the  afterbody  at  later  times.  It  may  be  useful 
to  admit  a  refinement  of  the  rear  complement  of  cells  to  retain  resolution  of 
the  chamber  geometry  as  well  as  to  reduce  this  disparity  in  cell  size. 

The  pressure  contours  show  that  in  spite  of  the  non-uniform  geometry  and 
the  flow  field,  the  distributions  quickly  become  nearly  one -dimensional  except 
right  near  the  base  of  the  afterbody.  Some  small  pressure  wiggles  are  apparent 
in  that  region.  It  is  not  presently  known  to  what  extent  these  reflect  the 
limited  complement  of  cells  or  the  local  cloistering  of  propellant  particles  due 
to  the  Helmholtz  mechanism  discussed  in  the  previous  chapter. 
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Figure  4.5  Porosity  Distributions  in  IZOsnui  Tank  Gun  at  Four  tiiiies 
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Figure  4.6  Ges-Pha.ee  Plow  Field  in  120inm  Tank  Gun  at  1.5  ms. 
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Figure  4.7  Ga^-PhsiSE  Flow  Field  in  X20Mn  Tank  Gun  at  3p0  msi 
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Figure  4,8  Gas-Phase  Flow  Field  in  120mm  Tank  Gun  at  4.5  ms. 
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Figure  4.9  Pressure  CDntcurs  in  120jnm  Tank  Gun  at  1.5  mSi 
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Figure  4.10  Pressure  Contours  in  IZOtnm  Tank  Gun  at  3.0  ms. 
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Figure  4.11  Pus  sure  Contours  in  120inm  Tank  Gun  at  4.5  ms. 


5.0  CONCLUSIONS 


We  have  developed  and  described  algorithms  to  support  two  key  modules  of 
the  NGEN  Code.  We  have  used  LCPFCT  as  the  basis  for  the  continuum  flow  solver 
(CFS)  and  developed  a  Large  Particle  Integrator  (LPI)  to  track  the  propellant 
grains  In  solid  propellant  charges. 

The  CFS  uses  LCPFCT  In  a  tlmespllt  formulation  which  Is  easy  to 
understand  and  modify.  We  have  presented  a  method  of  accommodating  non-uniform 
geometry  In  a  structured  mesh,  using  a  characteristic  based  analysis.  This 
analysis  permits  the  tlmestep  to  be  based  on  the  Courant  stability  condition 
for  Interior  cells  even  when  the  boundary  cells  will  violate  It  by  two  orders 
of  magnitude .  The  characteristic  based  analysis  has  been  formulated  In  a 
fashion  which  Is  simple  to  Implement  and  which  Imposes  a  minimal  computational 
burden.  The  approach  taken  to  the  treatment  of  non-uniform  geometry  Is  thought 
to  admit  massive  parallelization,  although  the  special  treatment  of  boundary 
cells  suggests  that  optimal  computer  architecture  may  be  such  that  Individual 
processors  Integrate  a  string  of  cells  rather  than  just  one.  More  Importantly, 
the  method  Is  thought  to  address  the  complexities  associated  with  the 
representation  of  Internal  boundaries  defined  by  Increment  containers.  The  LPI 
Is  also  simple  In  structure  and  provides  direct  control  of  the  porosity 
distribution  at  the  boundaries  of  charge  Increments. 

We  have  Illustrated  operability  of  the  algorithm  by  reference  to  two 
solid  propellant  data  bases,  one  for  a  155  mm  howitzer  and  the  other  for  a 
120  mm  tank  gun.  Each  calculation  required  about  10  minutes  on  an  SGI  Indigo 
work  station,  showing  the  algorithm  to  be  reasonably  fast  even  though  no 
special  measures  have  yet  been  taken  to  optimize  run  times. 


59 


REFERENCES 


Boris,  J.P.  and  Book,  D.L.  "Solution  of  Continuity  Equations  by  the 
Method  of  Flux-Corrected  Transport" 

Methods  in  Computational  Physics  pp  85-129  1976 

Boris',  J.P.,  Landsberg,  A.M. ,  Oran,  E.S.  and  Gardner,  J.H. 

"LCPFCT  -  A  Flux- Corrected  Transport  Algorithm  for  Solving  Generalized 
Continuity  Equations" 

NRL  Memorandxun  Report  NRL/MR/6410-93-7192,  Naval  Research  Laboratory  1993 

Chen,  J.L. ,  Kuo,  K.K,  and  Cheung,  F.B.  "Theoretical  Modeling  of  the 
Interior  Ballistic  Processes  in  an  Electrothermal  Gun" 

Proceedings  of  the  27th  Jannaf  Combustion  Meeting  1990 

Coffee,  T.P.  "A  Two-Dimensional  Model  for  the  Combustion  Chamber/ 

Gun  Tube  of  a  Concept  VI  C  Regenerative  Liquid  Propellant  Gun" 

Proceedings  of  the  27th  Jannaf  Combustion  Meeting  1990 

Cook,  D.C.,  Dyvlk,  J.A.  and  Chryssomallls ,  G.S.  "A  Multidimensional 
Electrothermal  Model"  Proceedings  of  the  26th  Jannaf  Combustion  Meeting  1989 

Courant,  R  and  Hilbert,  D.  "Methods  of  Mathematical  Physics" 

Interscience  1953 

Fisher,  E.B.  and  Graves,  K.W.  "Mathematical  Model  of  Double  Base 
Propellant  and  Ignition  in  the  81mm  Mortar"  Calspan  Report  DG-3029-D-1  1972 


Fitt,  A.D. ,  Crowley,  A.B.,  Aston,  J.A.G.  and  Toro,  E.F. 

"Contrasting  Numerical  Methods  for  Two-Dimensional  Two-Phase 
Interior  Ballistic  Flows" 

Proceedings  of  the  11th  International  Symposium  on  Ballistics  1989 

Gelperin,  N.l.  and  Einstein,  V.G.  "Heat  Transfer  in  Fluidized  Beds." 

Fluidization,  edited  by  Davidson,  J.F.  and  Harrison,  D.  Academic 

Press,  NY  1971 

Gibeling,  H.J.  and  McDonald,  H.  "An  Implicit  Numerical  Analysis 
for  Two-Dimensional  Turbulent  Interior  Ballistic  Flows" 

Ballistic  Research  Laboratory  Contract  Report  ARBRL-CR- 00523  1984 

Gough,  P.S.  "Modeling  of  Rigidized  Gun  Propelling  Charges" 

Ballistic  Research  Laboratory  Contract  Report  ARBRL-CR-00518  1983 

Gough,  P.S.  "The  XNOVAKTC  Code" 

Paul  Gough  Associates  Report  PGA-TR-86-1  1986 

Gough,  P.S.  "Formulation  of  a  Next-Generation  Interior  Ballistic 

Code"  Proceedings  of  the  29th  Jannaf  Combustion  Meeting  1991 


61 


Gough,  P.S.  "Formulation  of  a  Next-Generation  Interior  Ballistic 
Code  I.  Evaluation  of  Computational  Algorithms" 

Proceedings  of  the  30th  Jannaf  Combustion  Meeting  1992 

Groenenboom,  P.H.L.  and  Thomsen,  P.  "Two  Phase  Flow  Computation  for 
a  40mm  Granular  Propellant  Gun  in  Comparison  with  Experimental  Results" 
Proceedings  of.  the  11th  International  Symposium  on  Ballistics  1989 

Hadley,  G.  "Linear  Algebra"  Addlson-Vesley  1961 

Hsiao,  C.C.,  Phillips,  G.  and  Su,  F.  "First  Principles  Modeling  of  a  DNA 
60mm  ETC  Gun  Design"  Jannaf  Workshop  on  ETC  Modeling  and  Diagnostics  1991 

Kashiwa,  B.A. ,  Padlal,  N.T.  and  Butler,  D. 

"Toward  a  Comprehensive  Model  for  Electrothermal  Gun  Performance" 

Proceedings  of  the  27th  Jannaf  Combustion  Meeting  1990 

Krier,  H. ,  Shimpi,  S.A.  and  Adams,  M.J.  "Interior  Ballistic 
Predictions  Using  Data  From  Closed  and  Variable  Volume  Simulators." 

Univ.  of  Illinois  at  Urbana- Champaign.  TR-AAE-73-6  1973 

Landsberg,  A.M. ,  Boris,  J.P.,  Young,  T.R.  and  Scott,  R.J. 

"Computing  Complex  Shocked  Flows  Through  the  Euler  Equations" 

Naval  Research  Laboratory  1993 

Lohner,  R.  "Adaptive  H-refinement  on  3-D  Unstructured  Grids  for 

Transient  Problems"  AIAA  Paper  89-0653  1989 

Love,  A.E.A.  and  Pidduck,  F.B.  "Lagrange  Ballistic  Problem" 

Phil.  Trans.  Roy.  Soc.  vol  222  pp  167-226  1921 

May,  I.W.  and  Horst,  A.W.  "Charge  Design  and  Pressure  Waves  in 
Guns"  Papers  in  Astronautics  and  Aeronautics,  vol. 68,  Interior 
Ballistics  of  Guns,  Edited  by  H.  Krier  and  M.  Summerfield.  1979 

Meineke ,  E,  and  Heiser,  R.  "A  Complete  Numerical  Solution  of  the 
Interior  Ballistics  Chambrage  Problem" 

Proceedings  of  the  11th  International  Symposium  on  Ballistics  1989 

Oran,  E.S.,  Boris,  J.P.  and  Brown,  E.F.  "Fluid-Dynamic  Computations 
on  a  Connection  Machine  --  Preliminary  Timings  and  Complex  Boundary 
Conditions"  AlAA  Paper  90-0335  1990 

Robbins,  F.W.  "Comparison  of  TDNOVA  Results  with  an  Analytic  Solution" 
ARBRL-MR- 03299,  USA  ARRADCOM,  Ballistic  Research  Laboratory  1983 

Schmitt,  J.A.  and  Mann,  T.L. 

"An  Evaluation  of  the  Alpha  Code  in  its  One-Phase  Mode" 

ARBRL-MR-03081 ,  USA  ARRADCOM,  Ballistic  Research  Laboratory  1981 


62 


Schmitt,  J.A.  "A  Numerical  Algorithm  for  the  Multidimensional, 

Multiphase,  Viscous  Equations  of  Interior  Ballistics"  Proceedings 

of  the  2nd  Army  Conference  on  Applied  Mathematics  and  Computing  1984 

Slnha,  N. ,  Hosangadl,  A.  and  Dash,  S.M.  "Development  of  an  Upwind/ 

Implicit  Computational  Model  for  the  Advancement  of  Army  ETC  Guns" 

Jannaf  Workshop  on  ETC  Modeling  and  Diagnostics  1991 

Steffens,  U. ,  Krassln,  D.  and  Rlttel,  R.  "Two-Phase  Flow  Predictions  of 
Regenerative  Liquid  Propulsion  -  A  Comparison  with  Experimental  Results" 
Proceedings  of  the  Tenth  International  Symposium  on  Ballistics  1987 

Steffens,  U. ,  Rlttel,  R. ,  Lenselink,  H.  and  Florle,  K.  "Applications  of 
a  Three-Dimensional  Gasdynamlc  Simulation  Model  to  Interior  Ballistics" 
Proceedings  of  the  Eleventh  International  Symposium  on  Ballistics  1989 

Thomas,  P.D.  and  Lombard,  C.K.  "Geometric  Conservation  Law  and  its 
Application  to  Flow  Computations  on  Moving  Grids" 

AIAA  Journal  n.l7  pp  1030-1036  1979 

Thompson,  J.F.,  Varsl,  Z.U.A.  and  Mas ten,  C.V. 

"Numerical  Grid  Generation"  North  Holland,  New  York  1985 

Williams,  F.A.  "Combustion  Theory"  Addison-Wesley  1965 

Wlnsor,  N.K.  and  Goldstein,  S.A.  "Finite  Element  Modeling  of  an 
Electrothermal  Gun"  Proceedings  of  the  27th  Jannaf  Combustion  Meeting  1990 


63 


Intentionally  left  blank. 


64 


NOMENCIATDRE 


A  Area  of  computational  cell  in  LCPFCT 
a  Rate  of  propagation  of  Intergranular  disturbances 

fi  Bum  rate  coefficient 

b  Covolume 

c  Isentroplc  speed  of  sound  in  continuous  phase 
Cv  Specific  heat  at  constant  volume 

Cp  Specific  heat  at  constant  pressure 

Initial  diameter  of  grain 
Dp^  Effective  diameter  of  grain 

d  Surface  regression 

do  Initial  perforation  diameter  of  grain 

E  Total  energy  per  unit  volume  of  continuous  phase 

e  Internal  energy  of  continuous  phase 

e^g  Chemical  energy  of  propellant 

ep  Chemical  energy  of  igniter 

Interphase  drag  due  to  i-th  discrete  phase 

fj.  Radial  component  of  interphase  drag 
fg  Friction  factor 

fg  Axial  component  of  interphase  drag 

go  Constant  used  to  reconcile  units 

H  Thermal  parameter  for  cubic  profile  analysis  of  propellant  temperature 
h  Film  coefficient 

k  The  rmal  c onduc  t ivi ty 

L  Length  of  propellant  grain 

m  Mass  generation  per  unit  volume  due  to  propellant  combustion 

mij  Mass  generation  per  unit  volume  due  to  igniter  combustion 

N  Number  of  perforations  in  propellant  grain 
Ng  Number  of  continuous  phases 

N^  Number  of  dispersed  phases 

NUp  Nusselt  number 

n  Burn  rate  exponent 
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Number  density  of  i-th  dispersed  phase 

Pr  Prandtl  number 

p  Pressure 

q  Heat  transfer  per  unit  surface  area 

Heat  transfer  per  unit  volume  to  i-th  dispersed  phase 

qp  Heat  transfer  per  unit  volume  to  propellant 

Rep  Reynolds  number 

r  Radial  coordinate 

Surface  area  of  member  of  i-th  dispersed  phase 

Sp  Surface  area  of  a  propellant  grain 
T  Temperature  of  continuous  phase 

Temperature  of  i-th  dispersed  phase 

Tp  Surface  temperature  of  propellant 

t  Time 

u  Velocity  of  continuous  phase 

Velocity  of  i-th  dispersed  phase 
u  Axial  velocity  component  of  continuous  phase 

Up  Axial  velocity  component  of  propellant 

Volume  of  a  member  of  i-th  dispersed  phase 

Vp  Volume  of  a  propellant  grain 

Vj^  Diffusion  velocity  of  i-th  component  of  continuous  phase 
V  Radial  velocity  component  of  continuous  phase 

Vp  Radial  velocity  component  of  propellant 

Mass  fraction  of  i-th  component  of  continuous  phase 
z  Axial  coordinate 

Greek  Symbols 

a  Porosity 

Qj  Volume  fraction  of  i-th  component  of  continuous  phase 
Volume  fraction  of  i-th  component  of  dispersed  phase 

7  Ratio  of  specific  heats 

A  Volume  of  computational  cell  in  LCPFCT 
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(i  Viscosity 

p  Density  of  continuous  phase 

Pp  Density  of  propellant  grain 

Stress  tensor  for  i>th  dispersed  phase 
o  Intergranular  stress 

r  Shear  stress  tensor 

Rate  of  production  of  i>th  component  of  continuous  phase  due  to 

chemical  reactions 
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APPENDIX  A: 

CHARACTERISTIC  ANALYSIS  OF  COMTNUM  EQUATIONS 
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Our  interest  in  this  topic  is  principally  motivated  by  the  numerical 
ramifications  of  the  theory  of  characteristic  surfaces.  In  the  case  of  one 
dimensional  unsteady  flow,  the  existence  of  real  characteristic  directions 
enables  one  to  replace  the  system  of  partial  differential  equations  by  an 
equivalent  system  of  ordinary  differential  equations  in  which  the  derivatives 
are  taken  along  the  characteristic  lines.  When  we  proceed  to  a  larger  number 
of  independent  variables  an  analogous  result  holds  for  hyperbolic  systems  of 
equations.  Given  n  Independent  variables,  a  hyperbolic  system  is  one  that 
admits  the  existence  of  a  hypersurface  -of  dimension  n  -  1  such  that  only 
derivatives  interior  to  the  surface  appear  in  the  equations. 


We  proceed  as  follows.  In  Section  A.l,  we  discuss  the  theory  in  general 
for  a  quasilinear  system  of  partial  differential  equations  which  depend  on 
three  independent  coordinates.  Then,  in  Section  A. 2  we  deduce  the 
characteristic  forms  for  an  inviscid  two  dimensional  single  phase  flow  with 
non-homogeneous  source  terms.  By  treating  the  additional  terms  associated  with 
diffusion  as  embedded  in  the  non*homogeneous  source  terms  we  can  easily 
establish  "pseudocharacteristic"  forms,  which  are  of  value  in  the  formulation 
of  numerical  solution  algorithms  for  the  Navier-Stokes  equations.  Similarly, 
with  suitable  interpretation,  we  may  apply  the  results  obtained  here  in  the 
context  of  multiphase  flows. 

A.l  General  Formulation  of  Characteristic  Analysis 


Consider  a  system  of  partial  differential  equations 

at  az  ar 


A. 1.1 


where  V  and  D  are  n-dimensional  colvunn  vectors  and  A,  B,  C  are  n  x  n  square 
matrices.  The  concept  of  a  characteristic  surface  follows  naturally  from  the 
consideration  of  an  initial  value  problem  posed  for  a  surface 
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^(t,  z,  r)  -  0 


or,  in  general,  the  family  of  svirfaces  generated  by  the  parameter  4>o  such  that 


d(t,  z,  r)  -  0. 


A. 1.2 


Let  a,  he  coordinates  internal  to  the  surface  ^  bhen  ^  Itself  serves  as 
a  normal  coordinate.  On  4>  4>o  assume  that  ve  are  given  values  of  ^  and 
hence ,  values  of  and  .  The  surface  ^  is  said  to  be  free  if  Eqxiation 
A.  1.1  permits  the  determination  of  the  normal  derivative  and  characteristic 
if  it  does  not.  If  ^  is  characteristic  it  follows  that  A.  1.1  may  be 
expressed  in  terms  of  derivatives  with  respect  to  a  and  ^  alone,  that  is  to 
say,  derivatives  internal  to  the  characteristic  surface. 

By  means  of  the  chain  rule  for  differentiation,  A.  1.1  may  be  transformed 
into  derivatives  with  respect  to  a  and 


[Adt  +  «  D  -  [AOt  ♦  Bq,  +  COjlVa  -  (A^t  ♦ 


A. 1.3 


Accordingly,  the  question  of  whether  ^  is  free  or  characteristic  is  settled  by 
the  rank  of  the  matrix 


A  =  Adt  +  B^j  +  Cdp 


A. 1.4 


If  Rank  (A)  ^n,  the  system  A.  1.3  always  has  a  unique  solution  and  the 
surface  d  "  do  is  free.  However,  if  the  value  of  Rank  (A)<n  then  A.  1.3  does 
not  possess  a  solution  for  arbitrary  initial  data  on  the  surface  d  ”  do •  in  the 
latter  case  the  partial  differential  equation  represents  a  constraint  on  the 
data  as  expressed  by  the  condition  of  solvability  of  A. 1.3.  Thus  if  we  let  A* 
be  the  augmented  matrix  formed  by  appending  to  A  the  column  vector 
corresponding  to  the  right  hand  side  of  A. 1.3 

A*  =  [A;  D  -  [Aa^  +  Ba,  ♦  Cajd,  - 
Then  the  condition  of  solvability  is  (Hadley,  1961) 


Rank  (  A* )  ■  Rank  ( A ) 


A. 1.5 
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The  condition  Rank(d}<n  will  lead  to  a  partial  differential  equation 
for  ^(t,z,r)  In  the  form 

F(t,z,r,\&,^t,.^*.^r)  “  0 

It  being  assumed  that  A,B,C  are  functions  only  of  t,z,r  and-^. 
element  of  A  Is  a  homogeneous  linear  combination  of  4>t>4>z  and  4>x  it 
F  Is  homogeneous  of  order  k  >  1  in  these  quantities  so  that 

=  A*F(t,  z.  r. 

Accordingly  It  follows  that 

Because  of  the  degree  of  freedom  Induced  by  the  homogeneity  of  F  It  Is 
convenient  In  many  cases  to  append  an  additional  condition  corresponding  to  the 
normalization  of  the  vector  practice  the  most  convenient  choice 

Is  to  set  =  -1 .  This  corresponds  to  having  ^  In  the  form 


A. 1.6 

Since  every 
follows  that 

A. 1.7 

A. 1.8 


^  t(z,  r)  -  t 


A. 1.9 


so  that  =  -  1,  4>z  =  dt/dz  and  4>x  =  flt/3r. 


It  Is  useful  to  Interpret  these  results  geometrically  (Courant  and 
Hilbert,  1953).  We  may  think  of  4>  =  4>o  defining  a  surface  with  normal  vector 
proportional  to  Then  A.  1.6,  the  partial  differential  equation  for 
the  characteristic  surface.  Imposes  an  algebraic  constraint  on  the  components 
of  the  normal  at  each  point  In  the  (t,z,r)  space.  At  each  point  A. 1.6  defines 
a  family  of  planes  such  that  the  characteristic  surface  must  be  tangent  to  one 
of  them.  If  F  Is  not  linear,  the  envelope  of  this  family  of  planes  Is  a  cone, 
the  Honge  cone,  whose  generators  are  called  blcharacterlstlcs . 

According  to  A. 1.7  the  family  of  allowable  normal  vectors  at  a  given 
point  lies  on  the  surface  of  a  cone  whose  apex  is  the  point  in  question.  Thus 
A. 1.8  asserts  that  the  vector  is  perpendicular  not  only  to  the 
surface  of  normal  vectors,  but  also  to  the  vectors  themselves  since  they  lie 
along  the  cone.  In  fact,  the  vector  (F^^,  defines  a  blcharacterlstlc 
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direction.  From  A. 1.8  it  is  evident  that  it  must  lie  in  a  tangent  plane  of  a 
characteristic  surface.  However,  the  bicharacteristic  may  be  thought  of  as  the 
limit  of  the  line  of  intersection  of  neighboring  tangent  planes.  Thus  if  we 
write  b  as  a  bicharacteristic,  n  as  the  normal  to  a  tangent  plane  and  A  as  a 
parameter  which  labels  the  planes  at  a  given  point,  it  'follows  that 

b  B  n  X  [n  ^dAl.  Therefore  b  »  n  x  dA .  Since  both  n  and  lie  on  a 
I  dA  J  dA  dA 

tangent  plane  of  the  cone  defined  by  A.  1.6  we  see  that  b  is  parallel  to 
(F^,F^^,F^r)  •  bicharacteristic  ray  may  be  written  as 


A. 1.10 


This  result  may  be  used  to  eliminate  and  from  A.  1.6  and  to 

describe  the  characteristic  surface  by  reference  to  the  bicharacteristics. 

A. 2  Two-Dimensional  Inviscid  Single-Phase  Flow 

We  now  wish  to  deduce  complete  results  for  a  two-dimensional,  inviscid, 
single-phase  flow  with  a  local  source  term.  We  will  also  assume  that  the 
equations  are  subject  to  a  general  coordinate  transformation  for  computational 
purposes.  We  may  first  state  the  equations  in  cylindrical  coordinates  in  non¬ 
conservative  form  and  with  the  energy  equation  recast  to  eliminate  the  internal 
energy  in  favor  of  pressure  and  density  with  the  help  of  the  continuity 
equations  and  the  thermodynamic  identity 

frirl-ztl  ■ 

where  we  have  p,  pressure;  p,  density;  e,  internal  energy;  c,  isentropic  speed 
of  sound;  and  g^ ,  a  constant  to  reconcile  units. 


The  identification  of  the  non-homogeneous  terms  with  a  local  source  is 
for  illustrative  purposes  only.  As  noted  in  the  introduction,  the  formulism 
of  this  section  is  intended  to  be  applicable  to  more  complex  systems  of 
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equations  in  vhich  certain  partial  derivatives  may  also  be  treated  formally  as 
"non-homogeneous" .  Such  terms  will  Include  diffusion  terms  and  terms  related 
to  the  porosity  in  multiphase  flows. 


With  t,  time;  2,  axial  coordinate;  r,  radial  coordinate;  u,  axial 
velocity;  v,  radial  velocity;  and  m  a  local  source  yAiose  energy  content  is 
we  have 


2£  ♦pfiH 

Dt  [dz 


♦  izl  -  i  -  fZ  -  Cl 

8t)  r 


A. 2. 2 


P 


Du 

Dt 


-  mu 


A. 2. 3 


P 


Dt 


mv  =  {3 


A. 2. 4 


Dp  _  ^ 

Dt  g„  Dt 


m 


p(  —  )p 

3p 


♦ 


u*u 


e  -  p/p 


} 


u 


A. 2. 5 


and  where  D/Dt  is  the  convective  derivative. 

We  identify  these  equations  with  the  system 

=D 

at  az  ar 

by  setting 


A. 2. 6 
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u  p  0  0  - 

0  pu  0  go 

0  0  pu  0 

-  0  0  u 

So 


V  0  p  0 

0  pv  0  0 

0  0  pv  go 

-£!x  0  0  V 

So 


A. 2. 7 


At  this  point  we  have  established  the  balance  equations  in  a  form 
suitable  for  the  application  of  the  methodology  described  in  Section  A.l. 
However,  we  now  consider  a  transformation  of  coordinates  in  the  form 

r  =  t 

f  =  f  (t,  z,  r)  ,  A. 2. 8 

r)  =  Tj(t,z,r) 

We  assume  that  this  transformation  is  one  to  one  and  has  continuous  partial 
derivatives  and  that  a(f ,  fj)/3(z,  r)  0  so  that  we  can  also  write 

t  =  r 

z  =  z(T,f ,  fj)  ,  A. 2. 9 

r  =  r(T,  f.ij) 


We  have  used  a  separate  notation  r  for  the  time  in  the  transformed  frame.  This 
facilitates  the  use  of  subscripts  to  denote  partial  derivatives.  Thus  we 

understand  =  1—1  whereas  =  1—1  ^  is  an  arbitrary  property, 

latjz.r  larjf,. 
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Bearing  this  in  mind  we  define 
“  Zf  . 

A. 2. 10 

V,  - 

Thus  Ug,  and  are  the  velocity  components  In  a  cylindrical  coordinate  frame 
of  a  point  moving  so  that  It  Is  stationary  with  respect  to  the  transformed 
frame.  Evidently,  If  we  Impose  the  requirement  u„  ■  u  and  v,,  «  v  we  will  have 
selected  the  transformed  frame  to  coincide  with  a  Lagranglan  description  of  the 
fluid  whereas  the  choice  u,,  ■  v^  -  0  Implies  the  retention  of  ’  an  Eulerlan 
description,  possibly  In  a  different  coordinate  frame  established  by  a 
stationary  transformation. 


It  follows  that  the  balance  equations  A. 2. 6  subject  to  A. 2. 8  and  A. 2. 10 

become 

A  +  [  (B  -  AUg,)f,  +  (C  -  Avjfg]  +  [  (B  -  AuJfj,  +  (C  -  AvJrjJ  ^ 
dr  aj-  dr} 

A. 2. 11 

Thus  we  now  consider  the  characteristic  surfaces  for  the  system 

dr  ar  dT} 

where  we  Identify  and  C'  as: 

w  pr*  Pfr  0 

0  pw  0  gof, 

B' 

0  0  pw  g^fj 

-c^w/go  0  0  w 
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X 

PH, 

P^r 

0 

0 

px 

0 

go*?. 

0 

0 

px 

go*Jr 

*cVgo 

0 

0 

X 

and  where  we. have  Introduced 

w  =  (u  -  ♦  (v  -  , 

A. 2. 13 

X  «  (u  -  u„)»7,  +  (V  -  vjri^ 


Then  the  characteristic  surfaces  are  such  that  the  rank  of  A  is 

less  than  four  where 

A  ■  k4>f  *  .  A ,2. 14 

Thus  we  have 


A 


0 

0  0 

-c^^/go  0 


0  go(r*^f  +  Vz4>r,') 

P^  goCTr^f  ♦  Vt4>„) 

0  ^ 


and  where  we  have  introduced 


A. 2. 15 
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Now  in  our  applications  of  the  conditions  of  compatibility  we  shall 
always  require  that  either  »  0  or  that  ■  0  so  that  the  normal  lies  either 
in  the  r  -  f?  plane  or  the  r  -  f  plane.  We  assume  therefore  that  ■  0 .  The 
corresponding  results  for  the  case  ■  0  will  follow  from  considerations  of 
symmetry.  With  this  assumption  L  reduces  to 


A  ■ 


0 

0 


0  p'<t>  gofr^f 


A. 2. 16 


Now  we  consider  two  possibilities 

(i)  Let  ^  B  0.  Then  A  reduces  to 

0  P^fft  0 

0  0  0 

A  ■ 

0  0  0  g„r,^f 

0  0  0  0 

From  this  it  is  apparent  that  the  streamline 


^  ^  ®  0 

is  a  characteristic  direction.  In  order  to  establish  the  condition  of 
compatibility  we  now  introduce  a  and  /5  as  coordinates  internal  to  the 
characteristic  surface.  Then  we  may  write  the  augmented  matrix  as 
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i  -  (Opa  ♦  P[r,Of  ♦  ♦  P[fr“f  *  * 

:  -  {^P^  ♦  *  VtPv'^Vip  *  P[^A  * 

I 

:  ^2  -  <  P“U«  ♦  go[r*af  ♦  »7*0jPa  > 

:  -  {p0Up  *  Boi^tPc  *  »7*^,]Pp* 

:  ?3  -  t  pav„  ■*■  golfr^f  *  ^ 

•  -  {p/3Vp  +  go[rr/3f  +  Vj^]V0^ 

5  ?*  -  o[Po.  -  cVa/Bo  ]  -  PiVe-  C^Pp/go  ] 


where  a  and  are  defined  by  analogy  with  ^  . 


Accordingly,  if  we  write  ^^,1^1 . ,4  to  denote  the  members  of  the 

fifth  column  of  £^* ,  the  conditions  of  solvability  yield  the  following 
conditions  of  compatibility. 

r»^3  -  rr?2  =  0  .  A. 2. 17 


el 


A. 2. 18 


A  convenient  choice  for  a 


Q  «  r 


P  is 
^ « j? 


A. 2. 19 


and  we  will  adhere  to  this  convention.  The  use  of  a  separate  nomenclature  for 
the  coordinates  internal  to  the  characteristic  surface  is  again  motivated  by 
the  desire  to  maintain  clarity  in  respect  to  the  representation  of  the  partial 
derivatives . 
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Evidently 


•  a„  -  0 

•  “f  •  ^  -  1 

and  therefore, 

a  »  w  , 


^  “  X 


Then,  using  the  chain  rule  for  differentiation  and  noting  «  0  we  can  show  for 
any  variable  s  that 


Since  ^  -  0  implies  4>{/4>t  ~  Equation  A. 2. 17  may  be  identified  as  a  linear 

combination  of  the  two  momentum  equations.  Furthermore,  Equation  A.. 18  becomes 

Pt  *  wPf  -  —  (Pt  *  wPf )  =  “  X  fpf,  -  —  Pi,]  •  A. 2. 21 

go  I  go  J 


This  is  now  recognized  as  the  familiar  one  dimensional  result  with  the 
q -derivatives  taken  to  the  right  hand  side  and  treated  formally  as  non- 
homogeneous  terms.  Strictly  speaking  we  should  consider  the  singular  case  when 
w  B  0  so  that  ^  0  and  the  choice  A. 2. 19  cannot  be  made.  However,  the  same 
result  is  obtained  and,  in  any  case,  is  not  of  interest  here. 


(11)  Now  let  ^  f*  0 .  Then  perform  successively  the  following  row  and  column 
operations  to  A  as  given  by  Equation  A. 2. 16.  Add  gg/c^  times  column  one  to 
column  four;  subtract  times  row  two  from  row  one;  subtract  go^j^f/P^ 
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times  column  two  from  column  four;  subtract  times  column  three  from 

column  four;  subtract  times  row  three  from  row  one;  add  g^/c^  times  row 

four  to  row  one.  Then  A  is  equivalent  to  A  where 


Accordingly,  ^  can  only  be  characteristic  if 


That  is 


A. 2. 22 


W  ±  C 


Hz  *  fr  “  0 


A. 2. 23 


Then  identifying  Equation  A. 2. 22  with  F  in  Equation  A.  1.10  we  see  that  the 
bicharacteristics  satisfy  the  familiar  one  dimensional  form 

<*'  ■  - ^ — 7 -  A. 2. 24 

The  corresponding  condition  of  compatibility  is  easily  seen  to  be 


^  ^  So  ^  Q 

4>  ^  ^  C2 


Then  choosing  a  and  P  as  before,  Equation  A. 2. 19,  and  observing  df/d  ~  ^  1/c. 
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where  c,  «  c(ff  +  the  condition  of  compatibility  may  be  expressed  as 


/fF^TT 


(xp^  +  pfj.u^  ♦  P»?fV^)] 


(PXWp  ♦  goPp»?*)  ] 


±  [C3  “  (PXVp  +  goPp»?r)  ] 


A. 2. 25 


As  before,  the  ^-derivatives  are,  in  effect,  t;- derivatives  which  is  to  say 
derivatives  along  a  coordinate  curve  which  we  may  align  with  a  computational 
boundary.  Again,  Equation  A. 2. 25  is  analogous  with  the  one  dimensional  result 
with  the  cross  derivatives  (^-derivatives)  treated  formally  as  non-homogeneous 
terms . 


In  the  present  report  we  will  require  the  result  for  a  radial  boundary. 
Transforming  via  a  ^  r)  and  w  -*  x  we  have 


Pp  ±  ['^z'^0  *  ’Jr^p] 

goC. 


si 

go  U  ±  C,] 


A. 2. 26 


where  may  be  identified  from  Equation  A. 2. 25  as  the  non-homogeneous  group 
with  crossflow  terms  with  respect  to  a.  For  a  rectilinear  grid  we  may  restate 
Equation  A. 2. 26  as 


i£*Iv-v„±c)iP± 

at  ar 


pc 

go 


[  V  -  v„  ±  c  ) 


av 
ar . 


A. 2. 27 
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INTENTIONALLY  LEFT  BLANK. 
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APPENDIX  B: 

DESCRIPTION  OF  NGEN  INPUT  FILES 
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Intentionally  left  blank. 
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FILE  1  (10A8)  PROBLEM  TITLE  (1  CARD). 

TITLE  -  PROBLEM  TITLE,  UP  TO  80  ALPHANUMERIC  CHARACTERS. 

FILE  2  (515, 4X, nil)  LOGOUT  OPTIONS  (2  CARDS). 

NPRINT  -  IF  ZERO,  TABLES  OF  THE  STATE  VARIABLES  ARE  NOT 
PRINTED. 

-  IF  ONE,  TABLES  OF  THE  STATE  VARIABLES  ARE  PRINTED 
ON  A  LOGOUT  SCHEDULE  DETERMINED  BY  NSTEP  AND  DTLOG 
AS  DESCRIBED  IN  FILE  3. 

NSUMRY  -  IF  ZERO,  NO  SUMMARY  TABLES  ARE  PRODUCED  AT  THE 
CONCLUSION  OF  THE  RUN. 

-  IF  ONE,  SUMMARY  TABLES  ARE  PROVIDED  OF  THE  HISTORIES 
OF  THE  CONVENTIONAL  INTERIOR  BALLISTIC  VARIABLES 

NPLOT  -  IF  ZERO,  NO  ISOMETRIC  PLOTS  PRODUCED  ON  LOGOUT. 

-  IF  ONE,  PLOTS  OF  STATE  VARIABLES  PRODUCED  ON 
LOGOUT.  THESE  PLOTS  ARE  ISOMETRIC  VIEWS  OF  THE  STATE 
VARIABLES  AS  SELECTED  IN  ACCORDANCE  WITH  THE 
VALUES  OF  THE  ARRAY  IPLTV  DEFINED  BELOW.  PLOTTING  IS 
EFFECTED  BY  MEANS  OF  THE  POSTPROCESSOR  TGAP. 

NDSKW  -  IF  ZERO,  NO  DISC  STORAGE  ON  LOGOUT. 

-  IF  ONE,  SOLUTION  SAVED  ON  DISC  (UNIT  8)  ON  LOGOUT. 
NDSKR  -  IF  ZERO,  INITIAL  DISTRIBUTIONS  ARE  CONSTRUCTED  FROM 

INPUT  DATA. 

-  IF  NOT  ZERO,  INITIAL  DISTRIBUTIONS  ARE  READ  FROM 
UNIT  8  AND  CORRESPOND  TO  TIME  STEP  EQUAL  TO  NDSKR. 

IPLTV(I) ,1-1, . . . ,11  -  IF  IPLTV(I)-1,  THE  QUANTITY  TABULATED 
BELOW  WILL  BE  PLOTTED  AS  AN  ISOMETRIC  VIEW. 

OTHERWISE,  NOT. 

I  QUANTITY  PLOTTED  IF  IPLTV(I)-1. 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 


MESH. 

POROSITY. 

GRANULAR  STRESS. 

PRESSURE. 

DENSITY. 

GAS  AXIAL  VELOCITY. 

SOLID  AXIAL  VELOCITY. 

GAS  RADIAL  VELOCITY. 

SOLID  RADIAL  VELOCITY. 

GAS  TEMPERATURE. 

PARTICLE  SURFACE  TEMPERATURE. 
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FILE  3  (I5,F10.0,I5,F10.0,I5)  TERMINATION  AND  OUTPUT 

PARAMETERS  (1  CARD). 

NSTOP  -  INTEGRATION  STEP  FOR  TERMINATION. 

ZSTOP  -  PROJECTILE  DISPLACEMENT  FOR  TERMINATION  (CM) . 
NSTEP  -  NUMBER  OF  STEPS  BEFORE  PRINTOUT 
DTPRT  -  TIME  INTERVAL  FOR  PRINTOUT  (MSEC). 

NDEBUG  -  DEBUG  SWITCH.  IF  ZERO  THERE  IS  NO  DEBUG  PRINT. 
IF  NDEBUG  -  1,  DEBUG  PRINTING  OCCURS  ON  UNIT 
FORT. 50.  FILE  3.1  AND  3.2  ARE  REQUIRED. 

FILE  3.1  (315)  DEBUG  PARAMETERS  (1  CARD). 

NDTON  -  INTEGRATION  STEP  AT  WHICH  DEBUG  COMMENCES. 

NDTOFF  -  STEP  AT  WHICH  DEBUG  CEASES. 

NJLIS  -  NUMBER  OF  CELL  STRINGS. 

FILE  3.2  (315)  DEBUG  CELL  STRINGS  (NJLIS  CARDS). 

JDIR(I)  -  IF  1,  DEBUG  OCCURS  ON  AXIAL  SWEEP. 

IF  2,  DEBUG  OCCURS  ON  RADIAL  SWEEP. 

JLIS(I,1)  -  FIRST  RADIAL  OR  AXIAL  LOCATION  FOR 
DEBUG. 

JLIS(I,2)  -  LAST  RADIAL  OR  AXIAL  LOCATION  FOR 
DEBUG. 

DEBUG  PRINT  OCCURS  FOR  THE  LAST  5  AXIAL  OR  RADIAL 
CELLS  AT  EACH  RADIAL  OR  AXIAL  LOCATION  ACCORDING 
AS  JDIR  IS  EITHER  1  OR  2. 

FILE  4  (2F10. 0,315)  CONTINUUM  MESH  PARAMETERS  (1  CARD). 

CRN  -  COURANT  NUMBER,  MUST  BE  LESS  THAN  ONE.  A 
VALUE  OF  0.5  IS  RECOMMENDED. 

DIFFl  -  FRACTION  OF  ANTIDIFFUSIVE  FLUX  ALLOWED.  A 
VALUE  OF  0.999  IS  RECOMMENDED. 

NCELLR  -  TOTAL  NUMBER  OF  RADIAL  CELLS. 

NCHMIN  -  INITIAL  NUMBER  OF  AXIAL  CELLS  ALLOCATED  TO 
REGION  BEHIND  THE  BASE  OF  THE  AFTERBODY. 

NCHMAX  -  MAXIMUM  NUMBER  OF  AXIAL  CELLS  ALLOCATED  TO 
REGION  BEHIND  THE  BASE  OF  THE  AFTERBODY. 

THIS  DATUM  IS  PRESENTLY  INACTIVE. 

NABDYZ  -  NUMBER  OF  AXIAL  CELLS  ALLOCATED  TO  AFTERBODY. 

IF  >  0,  FILES  22  AND  23  ARE  REQUIRED. 

NABDYR  -  NUMBER  OF  RADIAL  CELLS  ALLOCATED  TO  BASE  OF 
AFTERBODY. 

NOTES:  (1)  NABDYZ  AND  NABDYR  ARE  ONLY  REQUIRED  IF  AN 
INTRUDING  AFTERBODY  IS  PRESENT. 

(2)  NCELLR  MUST  NOT  EXCEED  25. 

(3)  THE  SUM  OF  NCHMIN  AND  NABDYZ  MUST  NOT 
EXCEED  200. 
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FILE  5  (4F10.0)  AMBIENT  CONDITIONS  (1  CARD). 

PO  •  INITIAL  PRESSURE  OF  GAS  PHASE  (MPA) . 

TEMPO  -  INITIAL  TEMPERATURE  OF  GAS  PHASES  (DEG.K).  THIS  IS 
ALSO  THE  INITIAL  TEMPERATURE  OF  ANY  PROPELLANT  FOR 
WHICH  XTEMST  IS  ZERO  IN  FILE  8. 

GAMMAO  -  RATIO  OF  SPECIFIC  HEATS  OF  AMBIENT  GAS  (•). 

GMOLO  -  MOLECULAR  WEIGHT  OF  AMBIENT  GAS  (GM/GM-MOL) . 

BV  -  COVOLUME  (CM**3/GM) . 

FILE  6  (315)  CHARGE  REPRESENTATION  PARAMETERS  (1  CARD). 

NPRPS  -  NUMBER  OF  TYPES  OF  PROPELLANT  IN  CHARGE. 

MAXIMUM  OF  10. 

NBAGS  -  NUMBER  OF  INCREMENTS  OF  PROPELLANT.  MAXIMUM  OF 

LPIN,  A  PARAMETER  PRESENTLY  SET  EQUAL  TO  4.  LPIN 
SHOULD  NOT  BE  INCREASED-  TO  A  VALUE  LARGER  THAN  10. 
MODPRP  •  IF  0,  NO  SOLID  PROPELLANT  INCREMENTS  ARE  PRESENT. 
IF  1,  THE  PROPELLANT  IS  MODELED  ACCORDING  TO 
A  SIMPLE  LUMPED  PARAMETER  MODEL  IN  WHICH  THE 
CHARGE  CONSISTS  OF  A  SINGLE  INCREMENT  AND  EXPANDS 
ONLY  IN  THE  AXIAL  DIRECTION  AND  ACCORDING  TO  A 
LAGRANGEAN  (LINEAR)  DISTRIBUTION.  FILE  10.5  IS 
REQUIRED.  THE  PROPELLANT  IS  ASSUMED  TO  BE  TYPE  1. 

IF  2,  THE  INCREMENTS  ARE  MODELING  USING  A  TWO- 
DIMENSIONAL  MODEL.  FILE  11  IS  REQUIRED. 

NOTE  :  FILES  6.5,7,8,9,10  ARE  REPEATED,  AS  A  GROUP,  NPRPS  TIMES, 
ONCE  FOR  EACH  OF  THE  NPRPS  TYPES  OF  PROPELLANT  PRESENT 
IN  THE  CHARGE. 

IT  SHOULD  BE  NOTED  THAT  THE  SUBSCRIPT  USED  TO  DISTINGUISH 
THE  VARIOUS  TYPES  OF  PROPELLANT  IS  SUPPRESSED  IN  THE 
SUBSEQUENT  DISCUSSION. 

FILE  6.5  (10A8)  PROPELLANT  NAME  (1  CARD). 

PRNAME  -  PROPELLANT  DESCRIPTION.  UP  TO  80  ALPHANUMERIC 
CHARACTERS . 

FILE  7  (8F10.0)  SOLID  PHASE  CONSTITUTIVE  DATA  (1  CARD). 

XEE  -  SETTLING  POROSITY  OF  BED  (-). 

SEE  THE  NOTE  FOLLOWING  THE  DISCUSSION  OF  FILE  11 
CONCERNING  DEFAULT  PROPERTIES  OF  THIS  DATUM. 

XAPl  -  RATE  OF  PROPAGATION  OF  COMPRESSIVE  WAVE  IN  SETTLED 
BED  (CM/SEC). 

XAP2  -  RATE  OF  PROPAGATION  OF  UNLOADING  WAVE  (CM/SEC). 

XRHOP  -  DENSITY  OF  SOLID  PHASE  (GM/CC) . 

XKP  -  THERMAL  CONDUCTIVITY  OF  SOLID  PHASE  (J/CM-SEC-DEG .K) . 
XALFAP  -  THERMAL  DIFFUSIVITY  OF  SOLID  PHASE  (CM**2/SEC). 
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FILE  8  (4F10.0)  GAS  PHASE  CONSTITUTIVE  DATA  AND  INITIAL  PROPELLANT 
TEMPERATURE  (1  CARD). 

XGAM  -  RATIO  OF  SPECIFIC  HEATS  (-)• 

XGMOL  -  MOLECULAR  WEIGHT  (GM/GM-MOL) . 

XBV  -  COVOLUME  (CC/GM) . 

XTEMST  -  INITIAL  TEMPERATURE  OF  SOLID  PROPELLANT  (DEG.K).  IF 
XTEMST  IS  ENTERED  AS  ZERO,  IT  DEFAULTS  TO  TEMFO  AS 
DESCRIBED  IN  FILE  5. 

FILE  9  (I5/(8F10.0))  SOLID  PHASE  COMBUSTION  CHARACTERISTICS 

(2  OR  MORE  CARDS). 

NTB  '  NUMBER  OF  TABULAR  DATA  TO  DEFINE  BURN  RATE.  MAXIMUM 
OF  10. 

XTIG  -  IGNITION  TEMPERATURE  OF  SOLID  PHASE  (DEG.K). 

THIS  QUANTITY  STARTS  A  NEW  CARD. 

XECH  -  CHEMICAL  ENERGY  RELEASED  IN  COMBUSTION  (J/GM) . 

TMAXP(l)  -  MAXIMUM  PRESSURE  FOR  WHICH  CORRESPONDING 
COEFFICIENTS  ARE  APPLICABLE  IN  THE  LAW 
RD0T-TB1(1)+TB2(1)*P**TBN(1)  WHERE  P  IS 
PRESSURE  AND  RDOT  IS  REGRESSION  RATE. 

THIS  QUANTITY  STARTS  A  NEW  CARD. 

TBl(l)  -  BURN  RATE  ADDITIVE  CONSTANT  (CM/SEC). 

TB2(1)  -  BURN  RATE  PRE- EXPONENTIAL  FACTOR  (CM/SEC-MPA**BN) . 
TBN(l)  -  BURN  RATE  EXPONENT  (-). 


TMAXP(NTB)  -  MAXIMUM  PRESSURE  FOR  WHICH  CORRESPONDING 
COEFFICIENTS  ARE  APPLICABLE  IN  THE  LAW 
RD0T-TB1(NTB)+TB2(NTB)*P**TBN(NTB)  WHERE  P  IS 
PRESSURE  AND  RDOT  IS  REGRESSION  RATE. 

TBl(NTB)  -  BURN  RATE  ADDITIVE  CONSTANT  (CM/SEC). 

TB2(NTB)  -  BURN  RATE  PRE - EXPONENTIAL  FACTOR  ( CM/SEC -MPA**BN) . 
TBN(NTB)  -  BURN  RATE  EXPONENT  (-). 

NOTES:  (1)  A  NEW  CARD  IS  STARTED  FOR  TMAXP ( 1 ) , TMAXP ( 3 )  ETC,  BUT 
NOT  FOR  TMAXP(2) ,TMAXP(4)  ETC. 

(2)  IF  THE  PRESSURE  EXCEEDS  TMAXP (NTB) ,  THE  CORRESPONDING 
DATA  ARE  USED  AS  DEFAULT  VALUES. 

FILE  10  (15,6F10.0)  GRAIN  GEOMETRY  (1  CARD). 

NFORM  -  IF  ZERO,  GRAIN  IS  A  CYLINDER. 

-  IF  ONE,  GRAIN  IS  A  SPHERE. 

-  IF  TWO,  GRAIN  IS  AN  UNPERFORATED  STICK. 

XOD  -  EXTERNAL  DIAMETER  (CM) . 

XGLEN  -  LENGTH  (CM). 

XDPERF  -  DIAMETER  OF  PERFORATION  (CM). 

XNPERF  -  NUMBER  OF  PERFORATIONS  (-)• 
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FILE  10.5  (5F10.0)  DESCRIPTION  OF  INCREMENT  FOR  LUMPED  PARAMETER 

MODEL  WITH  LAGRANGEAN  MOTION  (1  CARD). 

NOTE:  FILE  REQUIRED  IF  AND  ONLY  IF  MODPRP  (FILE  6)  IS  EQUAL  TO  1. 

PROPMO  -  INITIAL  MASS  OF  PROPELLANT  (GM) . 

ZPROPl  -  POSITION  OF  FORWARD  BOUNDARY  OF  CYLINDRICAL  CORE  (CM) . 

ZPROP2  -  POSITION  OF  FORWARD  BOUNDARY  OF  PROPELLANT  (CM) .  THE 
POROSITY  IS  UNIFORM  FROM  ZERO  TO  ZPROPl  AND  TAPERS 
LINEARLY  TO  ONE  FROM  ZPROPl  TO  ZPROP2.  THESE  VALUES 
EXPAND  IN  PROPORTION  TO  THE  MOTION  OF  THE  PROJECTILE. 

NO  ACCOUNT  IS  TAKEN  OF  ANY  POSSIBLE  INTRUSION  BY  THE 
TUBE  OR  AFTERBODY  IN  CALCULATING  THE  POROSITY. 

RPROPl  -  POSITION  OF  OUTER  BOUNDARY  OF  CYLINDRICAL  CORE  (CM) . 

RPROP2  -  POSITION  OF  OUTER  BOUNDARY  OF  PROPELLANT  (CM) .  THESE 
DATA  CONTROL  THE  RADIAL  DISTRIBUTION  OF  POROSITY  IN 
A  FASHION  ANALOGOUS  TO  ZPROPl  AND  ZPROP2  EXCEPT  THAT 
THEY  REMAIN  CONSTANT  AT  ALL  TIMES. 

NOTE:  FILE  11  IS  REPEATED  NBAGS  TIMES,  ONCE  FOR  EACH  OF  THE  NBAGS 
INCREMENTS  OF  PROPELLANT. 

IT  SHOULD  BE  NOTED  THAT  THE  SUBSCRIPT  USED  TO  DISTINGUISH 
THE  VARIOUS  BAGS  OF  PROPELLANT  IS  SUPPRESSED. 

FILE  11  (3I5,5F10.0)  DESCRIPTION  OF  INCREMENT  FOR  TWO-DIMENSIONAL 

MODELING  (1  CARD). 

NOTE:  FILE  REQUIRED  IF  AND  ONLY  IF  MODPRP  (FILE  6)  IS  EQUAL  TO  2. 

MPRP  -  POINTER  TO  FILE  OF  PROPELLANT  PROPERTIES  WHICH 

CHARACTERIZE  THE  MAIN  CHARGE  CONTAINED  IN  THE  INCREMENT. 
MPRP  MUST  BE  GREATER  THAN  ZERO  AND  LESS  THAN  OR 
EQUAL  TO  NPRPS  (FILE  5). 

NZC  -  NUMBER  OF  AXIAL  PARTICLES  FOR  INCREMENT.  MUST  NOT 
EXCEED  PARAMETER  LPIZ,  PRESENTLY  EQUAL  TO  50. 

NRC  -  NUMBER  OF  RADIAL  PARTICLES  FOR  INCREMENT.  MUST  NOT 
EXCEED  LPIR,  PRESENTLY  EQUAL  TO  20. 

XCHWT  -  INITIAL  MASS  OF  MAIN  CHARGE  IN  INCREMENT  (GM) . 

XZR  -  INITIAL  POSITION  OF  REAR  BOUNDARY  (CM) . 

XZF  -  INITIAL  POSITION  OF  FORWARD  BOUNDARY  (CM) . 

XRI  -  INITIAL  POSITION  OF  INNER  BOUNDARY  (CM) . 

XRO  -  INITIAL  POSITION  OF  OUTER  BOUNDARY  (CM) . 

NOTES:  (1)  IF  XZR,  XZF,  XRI  OR  XRO  IS  ENTERED  AS  ZERO,  IT  IS 

DEFAULTED  TO  THE  CORRESPONDING  BOUNDARY  OF  THE  GUN 
CHAMBER. 

(2)  THE  INCREMENT  IS  TAKEN  TO  BE  UNIFORMLY  DISTRIBUTED 
WITHIN  THE  ENVELOPE  DEFINED  BY  THESE  DATA  LESS  ANY 
INTRUSIONS  BY  THE  TUBE  WALL  OR  THE  AFTERBODY. 

(3)  PRESENT  CODING  ASSUMES  MULTIPLE  INCREMENTS  TO  BE 
ARRANGED  END-TO-END. 
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FILE  12  (415)  IGNITER  DISCHARGE  TABLE  COUNTERS  AND  OPTIONS 

(1  CARD). 

NTABIG  -  IF  ZERO,  A  TABULAR  REPRESENTATION  OF  AN 

IGNITION  STIMULUS  VIEWED  AS  AN  EXTERNALLY 
INJECTED  SOURCE  IS  NOT  CONSIDERED. 

-  IF  ONE.  AN  EXTERNALLY  INJECTED  IGNITION 
SOURCE  IS  CONSIDERED.  VALUES  OF  JZP,  JRP 
AND  JTP  MUST  BE  SPECIFIED  AND  FILES  13,14,15 
16  AND  17  MUST  BE  INCLUDED. 

JZP  -  NUMBER  OF  AXIAL  STATIONS  IN  DISCHARGE  TABLE 
FOR  CASE  WHEN  NTABIG  EQUALS  ONE. 

JZP  MUST  NOT  EXCEED  TWENTY. 

JRP  -  NUMBER  OF  RADIAL  STATIONS  IN  DISCHARGE  TABLE 
FOR  CASE  WHEN  NTABIG  EQUALS  ONE. 

JRP  MUST  NOT  EXCEED  EIGHT. 

JTP  -  NUMBER  OF  TIME  LEVELS  IN  DISCHARGE  TABLE 
FOR  CASE  WHEN  NTABIG  EQUALS  ONE. 

JTP  MUST  NOT  EXCEED  TWENTY. 

FILE  13  (3F10.0)  ENERGY  OF  EXTERNAL  IGNITION  SOURCE  (1  CARD). 
NOTE:  THIS  FILE  IS  REQUIRED  IF  AND  ONLY  IF  NTABIG  IS 
EQUAL  TO  ONE. 

EIG  -  ENERGY  OF  IGNITER  GAS  (J/GM) . 

GAMIG  -  RATIO  OF  SPECIFIC  HEATS  OF  IGNITER  GAS  (-)• 

GMOLIG  -  MOLECULAR  WEIGHT  OF  IGNITER  GAS  (GM/GMOL) . 

FILE  14  (8F10.0)  AXIAL  POSITIONS  FOR  DISCHARGE  TABLE  (1  TO 

3  CARDS). 

NOTE:  THIS  FILE  IS  REQUIRED  IF  AND  ONLY  IF  NTABIG  IS  EQUAL 
TO  ONE. 

ZPHI(I) ,I-1,JZP  -  AXIAL  POSITIONS  (CM). 

FILE  15  (8F10.0)  RADIAL  POSITIONS  FOR  DISCHARGE  TABLE  (1  CARD). 
NOTE:  THIS  FILE  IS  REQUIRED  IF  AND  ONLY  IF  NTABIG  IS  EQUAL 
TO  ONE. 

RPHI(I),I-1,JRP  -  RADIAL  POSITIONS  (CM). 

FILE  16  (8F10.0)  TIME  LEVELS  FOR  DISCHARGE  TABLE  (1  TO  3  CARDS). 
NOTE:  THIS  FILE  IS  REQUIRED  IF  AND  ONLY  IF  NTABIG  IS  EQUAL 
TO  ONE. 

TPHI(I) ,I-1,JTP  -  TIME  LEVELS  (MSEC). 
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FILE  17  (8F10.0)  DISCHARGE  TABLE  (JRP*JTP  TO  3*JRP*JTP  CARDS). 
NOTE:  THIS  FILE  IS  REQUIRED  IF  AND  ONLY  IF  NTABIG  IS  EQUAL 
TO  ONE. 

PHI (1,1,1)  -  FIRST  VALUE  OF  RATE  OF  DISCHARGE  PER  UNIT  VOLUME 
(GM/CC-SEC). 

PHI (2, 1,1)  -  SECOND  VALUE. 


PHI(JZP,1,1) 
PHI(1,2,1)  - 


-  VALUE  AT  LAST  AXIAL  POSITION,  FIRST  RADIAL 
POSITION  AND  FIRST  TIME. 

VALUE  AT  FIRST  AXIAL,  SECOND  RADIAL  POSITION. 
THIS  ENTRY  STARTS  A  NEW  CARD. 


PHI(JZP,JRP,JTP)  -  LAST  VALUE. 

FILE  18  (2F10. 0,515)  PROJECTILE  POSITION,  MASS, 

AND  BORE  RESISTANCE  DATA  (1  CARD) . 

ZPO  -  INITIAL  LOCATION  OF  BASE  OF  PROJECTILE 
OR  SABOT  (CM).  THIS  CORRESPONDS  TO  THE 
FORWARD  BOUNDARY  OF  THE  COMPUTATIONAL 
DOMAIN  AND  THE  AFTERBODY,  IF  PRESENT,  IS 
VIEWED  AS  INTRUDING  INTO  THE  REGION  TO 
REAR  OF  THIS  LOCATION. 

PRMASS  -  PROJECTILE  MASS  (KG) . 

NBRES  -  NUMBER  OF  ENTRIES  IN  TABULAR  DESCRIPTION  OF 
BORE  RESISTANCE.  MUST  NOT  EXCEED  10. 

IBRES  -  TYPE  OF  LAW  FOR  BORE  RESISTANCE. 

IF  1,  RESISTANCE  GIVEN  DIRECTLY  BY  INTERPOLATION 
OF  TABULAR  DATA  OF  FILE  28. 

IF  2,  INTERPOLATED  VALUE  MULTIPLIED  BY 
7.2/V**0.6 

WHERE  V  IS  PROJECTILE  VELOCITY  IN  FT/SEC. 

IF  3,  INTERPOLATED  VALUE  MULTIPLIED  BY 
(1+.0004414V)/1+.005046V) 

WHERE  V  IS  PROJECTILE  VELOCITY  IN  IN/SEC. 

IF  11,12  OR  13,  THE  RESISTANCE  IS  AS  FOR  1,2  OR  3 
RESPECTIVELY  AND  THE  RESISTANCE  DUE  TO 
COMPRESSED  AIR  IN  FRONT  OF  THE  PROJECTILE 
IS  ALSO  TAKEN  INTO  ACCOUNT  ACCORDING  TO  AN 
ANALYTICAL  FORMULA  BASED  ON  STEADY-STATE 
SHOCK  THEORY. 
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FILE  19  (8F10.0)  BORE  RESISTANCE  TABLE  (1  TO  3  CARDS). 
ZBRES(l)  -  FIRST  VALUE  OF  PROJECTILE  DISPLACEMENT 
AT  WHICH  BORE  RESISTANCE  IS  SPECIFIED 
(CM). 

FBRES(l)  -  CORRESPONDING  VALUE  OF  BORE  RESISTANCE 
(MPA). 


ZBRES(NBRES)  -  LAST  VALUE  OF  DISPLACEMENT. 

FBRES(NBRES)  -  CORRESPONDING  VALUE  OF  BORE  RESISTANCE. 

FILE  20  (15)  TUBE  GEOMETRY  FILE  COUNTER  (1  CARD). 

NBYE  -  NUMBER  OF  PAIRS  OF  DATA  TO  DEFINE  TUBE  GEOMETRY. 
MINIMUM  OF  TWO  AND  MAXIMUM  OF  TEN. 

FILE  21  (2F10.0)  TUBE  GEOMETRY  DATA  (NBYE  CARDS). 

ZBYE(l)  -  FIRST  AXIAL  POSITION  (CM). 

RBYE(l)  -  CORRESPONDING  RADIAL  POSITION  (CM). 

ZByE(2)  -  SECOND  AXIAL  POSITION.  STARTS  A  NEW  CARD. 


RBYE(NBYE) 

NOTE;  THE  GEOMETRY  IS  ARBITRARY  EXCEPT  THAT  ZBYE(I)  MUST  BE 
GREATER  THAN  ZBYE(I*1)  AND  THAT  THE  RESULTING 
CONFIGURATION  MUST  NOT  CROSS  ANY  AXIAL  LINE  IN  MORE 
THAN  TWO  LOCATIONS.  PRESENT  CODING  IN  NGEN  ASSUMES 
THAT  THE  OPEN  CELLS  FOR  EACH  AXIAL  SWEEP  ARE  SINGLY 
CONNECTED.  CODING  ALSO  ASSUMES  THAT  THE  TUBE  WALL 
IS  MODERATELY  TAPERED.  TUBES  WITH  VERY  STRONG  CHAMBRAGE 
MAY  PRESENT  COMPUTATIONAL  DIFFICULTIES  WITH  THE  PRESENT 
ALGORITHM. 

FILE  22  (15)  AFTERBODY  FILE  COUNTER  (I  CARD). 

NOTE:  THIS  FILE  IS  ONLY  REQUIRED  IF  NABDYZ  (FILE  4)  IS  GREATER 
THAN  ZERO. 

NBYI  -  NUMBER  OF  PAIRS  OF  DATA  TO  DEFINE  AFTERBODY  GEOMETRY. 

MINIMUM  OF  TWO  AND  MAXIMUM  OF  TEN. 

FILE  23  (2F10.0)  AFTERBODY  GEOMETRY  (NBYI  CARDS). 

NOTE:  THIS  FILE  IS  ONLY  REQUIRED  IF  NABDYZ  (FILE  4)  IS  GREATER 
THAN  ZERO. 

ZBYI(l)  -  FIRST  AXIAL  POSITION  (CM). 

RBYI(l)  -  CORRESPONDING  RADIAL  POSITION  (CM). 

ZBYI(2)  -  SECOND  AXIAL  POSITION.  STARTS  A  NEW  CARD. 


RBYI(NBYI) 

NOTE:  THE  COMMENTS  CONCERNING  THE  TUBE  GEOMETRY  DATA  APPLY 
HERE  ALSO.  IN  ADDITION  NOTE  THAT  THE  AFTERBODY  DATA 
ARE  NOT  NECESSARILY  REFERENCED  TO  THE  BREECHFACE.  AN 
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INTERNAL  ADJUSTMENT  COLLOCATES  ZBYI(NBYI)  WITH  THE 
BASE  OF  THE  PROJECTILE/SABOT. 

ALSO  NOTE  THAT  IT  IS  ASSUMED  THAT  THE  BASE  OF  THE  AFTERBODY, 
LIKE  THAT  OF  THE  BASE  OF  THE  PROJECTILE/SABOT,  IS  FLAT  AND 
OF  SUFFICIENT  RADIAL  EXTENSION  TO  BE  COVERED  BY  AT  LEAST 
ONE  RADIAL  CELL  (SEE  NABDYR,  FILE  4).  FINALLY,  THE  EXPOSED 
PART  OF  THE  BASE  OF  THE  PROJECTILE/SABOT  MUST  BE  LARGE 
ENOUGH  TO  BE  COVERED  BY  AT  LEAST  3  RADIAL  CELLS  IN  THE 
PRESENT  VERSION  OF  THE  CODE. 
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Intentionally  left  blank. 
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NO.  OF 

COPIES  ORGANIZATION 

2  ADMINISTRATOR 
ATTN  DTIC  DDA 

DEFENSE  TECHNICAL  INFO  CTR 
CAMERON  STATION 
ALEXANDRIA  VA  22304-6145 

1  DIRECTOR 

ATTN  AMSRL  OP  SD  TA 
US  ARMY  RESEARCH  LAB 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1145 

3  DIRECTOR 

ATTN  AMSRL  OP  SD  TL 
US  ARMY  RESEARCH  LAB 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1145 

1  DIRECTOR 

ATTN  AMSRL  OP  SD  TP 
US  ARMY  RESEARCH  LAB 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1145 


ABERDEEN  PROVING  GROUND 

5  DIR  USARL 

ATTN  AMSRL  OP  AP  L  (305) 


97 


NO.  OF 

COPffiS  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


1  HQDA 

ATTN  SARD  TR  MS  K  KOMINOS 
PENTAGON 

W  ASHINGTON  DC  20310-0103 
1  HQDA 

ATTN  SARD  TR  DR  R  CHATT 
PENTAGON 

WASHINGTON  DC  20310-0103 

1  CHAIRMAN  DOD  EXPLOSIVES  SAFETY  BD 

HOFFMAN  BLDG  1  RM  856  C 
2461  EISENHOWER  AVE 
ALEXANDRIA  VA  22331-0600 

1  HQS  US  ARMY  MATERIEL  CMD 

ATTN  AMCICP  AD  M  FISETTE 
5001  EISENHOWER  AVE 
ALEXANDRIA  VA  22333-0001 

1  US  ARMY  BMDS  CMD 

ADVANCED  TECHLOY  CTR 
PO  BOX  1500 

HUNTSVILLE  AL  35807-3801 

1  OFC  OF  THE  PRODUCT  MGR 

ATTN  SFAE  AR  HIP  IP 
MR  R  DE  KLEINE 

155MM  HOWITZER  M109A6  PALADIN 
PCNTYARSNLNJ  07806-5000 

3  PM  ADV  FIELD  ARTLRY  SYSTEM 

ATTN  SFAE  ASM  AF  E 
LTC  A  ELLIS 
TKURIATA 
J  SHIELDS 

PCNTYARSNLNJ  07801-5000 

1  PM  ADV  FIELD  ARTLRY  SYSTEM 

ATTN  SFAE  ASM  AF  Q  W  WARREN 
PCNTYARSNLNJ  07801-5000 

1  CDR  US  ARMY  ARDEC 

ATTN  AMSMC  PBM  A  SIKLOSI 
PROD  BASE  MODERNIZATION  AGENCY 
PCNTYARSNLNJ  07806-5000 

1  CDR  US  ARMY  ARDEC 

ATTN  AMSMC  PBM  E  L  LAIBSON 
PROD  BASE  MODERNIZATION  AGENCY 
PCNTY  ARSNL  NJ  07806-5000 


1  PM  PEO  ARMAMENTS 
ATTN  AMCPMTMA 
TANK  MAIN  ARMAMENT  SYSTEM 
PCNTYARSNLNJ  07806-5000 

1  PM  PEO  ARMAMENTS 
ATTN  AMCPM  TMA  105 
TANK  MAIN  ARMAMENT  SYSTEM 
PCNTYARSNLNJ  07806-5000 

1  PM  PEO  ARMAMENTS 

ATTN  AMCPM  TMA  120 
TANK  MAIN  ARMAMENT  SYSTEM 
PCNTYARSNLNJ  07806-5000 

1  PM  PEO  ARMAMENTS 

ATTN  AMCPM  TMA  AS  H  YUEN 
TANK  MAIN  ARMAMENT  SYSTEM 
PCNTYARSNLNJ  07806-5000 

2  CDR  US  ARMY  ARDEC 
ATTN  SMCAR  CCH  V 
CMANDALA 
EFENNELL 

PCNTYARSNLNJ  07806-5000 

1  CDR  US  ARMY  ARDEC 

ATTN  SMCAR  CCH  T  L  ROSENDORF 
PCNTY  ARSNL  NJ  07806-5000 

1  CDR  US  ARMY  ARDEC 

ATTN  SMCAR  CCS 
PCNTYARSNLNJ  07806-5000 

1  CDR  US  ARMY  ARDEC 

ATTN  SMCAR  AEE  J  LANNON 
PCNTYARSNLNJ  07806-5000 

11  CDR  US  ARMY  ARDEC 

ATTN  SMCAR  AEE  B 
A  BEARDELL 
D  DOWNS 
S  EINSTEIN 
S  WESTLEY 
S  BERNSTEIN 
J  RUTKOWSKI 
B  BRODMAN 
P  O’REILLY 
R  CIRINCIONE 
P  HUI 
J  O’REILLY 

PCNTYARSNLNJ  07806-5000 
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NO.  OF 

COPIES  ORGANIZATibN 


NO.  OF 

COPIES  ORGANIZATION 


S  COMMANDER 

ATTN  SMCAR  AEE  WW 

MMEZGER 

J  PINTO 

D  WIBGAND 

PLU 

CHU 

US  ARMY  ARDEC 
PCNTYARSNLNJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  AES  S  KAPLOWTIZ 

US  ARMY  ARDEC 

PCNTY  ARSNL  NJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  HFM  E  BARRIERES 
US  ARMY  ARDEC 
PCNTYARSNLNJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  FSA  T  M  SALSBURY 
US  ARMY  ARDEC 
PCNTYARSNLNJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  FSA  F  LTC  R  RIDDLE 
US  ARMY  ARDEC 
PCNTYARSNLNJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  FSC  G  FERDINAND 
US  ARMY  ARDEC 
PCNTYARSNLNJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  FS  T  GORA 
US  ARMY  ARDEC 
PCNTYARSNLNJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  FS  DH  J  FENECK 

US  ARMY  ARDEC 

PCNTY  ARSNL  NJ  07806-5000 

3  COMMANDER 

ATTN  SMCAR  FSS  A 

RKOPMANN 

BMACHEK 

L  PINDER 

US  ARMY  ARDEC 

PCNTYARSNLNJ  07806-5000 


1  COMMANDER 

ATTN  SMCAR  FSN  N  K  CHUNG 
US  ARMY  ARDEC 
PCNTYARSNLNJ  07806-5000 

2  DIR  BENET  WEAPONS  LABS 
ATTN  SMCAR  CCB  RA 

G  P  O’HARA 
G  A  PFLEGL 

WATERVLIETNY  12189-4050 

1  DIR  BENET  WEAPONS  LABS 

ATTN  SMCAR  CCB  RT  S  SOPOR 
WATERVLIETNY  12189-4050 

1  DIR  BENET  WEAPONS  LABS 
ATTN  SMCAR  CCB  S  F  HEISER 
WATERVLIETNY  12189-4050 

2  CDR  US  ARMY  RSRCH  OFC 
ATTN  TECHNICAL  LIBRARY 
DMANN 

POBOX  12211 

RSCHTRIPKNC  27709-2211 

1  CDR  USACECOM 

ATTN  ASQNC  ELC  IS  L  R  MYER  CENTER 
R&D  TECHNICAL  LIBRARY 
FORT  MONMOUTH  NJ  07703-5301 

1  CMDT  US  ARMY  AVIATION  SCHOOL 

ATTN  AVIATION  AGENCY 
FORT  RUCKER  AL  36360 

1  PM  US  TANK  AUTOMOTIVE  CMD 

ATTN  AMCPM  ABMS  T  DEAN 
WARREN  MI  48092-2498 

1  PM  US  TANK  AUTOMOTIVE  CMD 

ATTN  SFAE  ASM  BV 
FIGHTING  VEHICLE  SYSTEMS 
WARREN  MI  48397-5000 

1  PM  ABRAMS  TANK  SYSTEM 

ATTN  SFAE  ASM  AB 
WARREN  MI  48397-5000 

1  DIR  HQ  TRAC  RPD 

ATTN  ATCD  MA 
FORT  MONROE  VA  23651-5143 


NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


1  COMMANDER 

ATTN  STRBE  WC 

US  ARMY  BELVOK  R&D  CTR 

FORT  BELVOIR  VA  22060-5006 

1  DIRECTOR 

ATTN  ATRC  L  MR  CAMERON 
US  ARMY  TRAC  FT  LEE 
FORT  LEE  VA  23801-6140 

1  COMMANDANT 

US  ARMY  CMD  &  GEN  STAFF  COLLEGE 
FORT  LEAVENWORTH  KS  66027 

1  COMMANDANT 

ATTN  REV  AND  TRNG  LIT  DIV 
US  ARMY  SPECIAL  WARFARE  SCHOOL 
FORT  BRAGG  NC  28307 

1  COMMANDER 

ATTN  SMCAR  QA  HI  LIB 

RADFORD  ARMY  AMMUNITION  PLANT 

RADFORD  VA  24141-0298 

1  COMMANDER 

ATTN  AMXST  MC  3 

US  ARMY  FRGN  SCIENCE  &  TECHLGY  CTR 
220  SEVENTH  STREET  NE 
CHRLTTESVLLE  VA  22901-5396 

1  COMMANDANT 

ATTN  ATSF  CD  COL  T  STRICKLIN 
US  ARMY  FIELD  ARTLRY  CTR  &  SCHOOL 
FT  SILL  OK  73503-5600 

1  COMMANDANT 

ATSF  CN  P  GROSS 

US  ARMY  FIELD  ARTLRY  CTR  &  SCHOOL 
FT  SILL  OK  73503-5600 

1  CMDT  US  ARMY  ARMOR  SCHOOL 
ATTN  ATZK  CD  MS  M  FALKOVITCH 
ARMOR  AGENCY 

FORT  KNOX  KY  40121-5215 

2  CDR  NAVAL  SEA  SYSTEMS  CMD 

ATTN  SEA  62R 

SEA  64 

WASH  DC  20362-5101 


1  CDR  NAVAL  AIR  SYSTEMS  CMD 

ATTN  AIR  954  TECH  LIBRARY 
WASH  DC  20360 

4  CDR  NAVAL  RSRCH  LAB 

ATTN  TECHNICAL  LIBRARY 
CODE  4410 
K  KAILASANATE 
J  BORIS 
EORAN 

WASH  DC  20375-5000 

1  OFFICE  OF  NAVAL  RSRCH 

ATTN  CODE  473  R  S  MILLER 
800  N  QUINCY  STREET 
ARLINGTON  VA  22217-9999 

1  OmCE  OF  NAVAL  TECHLGY 

ATTN  ONT  213  D  SIEGEL 
800  N  QUINCY  ST 
ARLINGTON  VA  22217-5000 

1  CDR  NAVAL  SURFACE  WARFARE  CTR 
ATTN  CODE  730 
SILVER  SPRING  MD  20903-5000 

1  CDR  NAVAL  SURFACE  WARFARE  CTR 

ATTN  CODE  R  13  R  BERNECKER 
SILVER  SPRING  MD  20903-5000 

7  CDR  NAVAL  SURFACE  WARFARE  CTR 

ATTN  T  C  SMITH 
KRICE 
S  MITCHELL 
S  PETERS 
J  CONSAGA 
C  GOTZMER 
TECHNICAL  LIBRARY 
INDIAN  HEAD  MD  20640-5000 

1  CDR  NAVAL  SURFACE  WARFARE  CTR 

ATTN  CODE  G30  GUNS  &  MUNITIONS  DIV 
DAHLGRENVA  22448-5000 

1  CDR  NAVAL  SURFACE  WARFARE  CTR 
ATTN  CODE  G32  GUNS  SYSTEMS  DIV 
DAHLGRENVA  22448-5000 

1  CDR  NAVAL  SURFACE  WARFARE  CTR 
ATTN  CODE  G33  T  DORAN 
DAHLGRENVA  22448-5000 
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